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SUMMARY 

A comparison of flow separation in transonic flows is made using various 
computational schemes which solve the Euler and the Navier-Stokes equations of 
fluid mechanics. The flows examined are computed using several simple two- 
dimensional configurations including a backward facing step and a bump in a 
channel. Comparison of the results obtained using shock fitting and flux 
vector splitting methods are presented and the results obtained using the 
Euler codes are compared to results on the same configurations using a code 
which solves the Navier-Stokes equations. 
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Chapter 1 


INTRODUCTION 

Although much progress has been made in the field of computational 
fluid dynamics (CFD), the computation of transonic flows containing 
shock waves and of flows exhibiting separation is still a challenge. 
The most commonly used methods to compute such flows solve either the 
time-dependent Euler or Navier-Stokes equations on a fixed grid. The 
Euler equations allow for rotational flows but neglect viscous effects, 
whereas the Navier-Stokes equations take into account viscous effects. 
Both sets of equations, when expressed in integral form, are correct 
even when discontinuities in the flow field are present. The 
compressible potential equation, which has been used extensively for 
aerodynamic prediction, does not allow for rotational flows and viscous 
effects and is therefore not considered in this study. However, when 
shock waves are weak and flow separation is not expected, the potential 
model may provide a good approximation. Furthermore, when coupled with 
boundary layer equations, viscous-lnviscid methods have been 
successfully used to compute flows with separation. 

The development of CFD and, more specifically, methods to solve 
compressible flow problems extends back into the 1950 ' s and a complete 
history is beyond the scope of this introduction. The mathematical 
theory of these numerical approximations has been developing rapidly as 
has been the mathematical theory of hyperbolic conservation laws and 
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shock waves. The 1948 book by Courant and Friedrichs [1]* and the 1973 
SIAM publication by Lax [2] were important contributions to the 
literature which provide much of the basic theory. The 1967 text by 
Richtmeyer and Morton [3] is now considered a classic text on the 
application of difference methods to initial value problems. More 
recently, the books by Smoller [4] and Majda [5] provide up-to-date 
exposition of the mathematical theory. While there have been many books 
on CFD published, the recent text by Anderson, Pletcher and Tannehill 
[6] is one of the most popular. Finally, the 1986 survey article by Roe 
contains an excellent discussion of the development and fairly recent 
state of affairs in numerical schemes for the Euler equations [7]. 

Early schemes were relatively simple explicit methods such as the 
Lax-Wendroff scheme [8] and its derivative, the two-step MacCormack 
Scheme [9]. These schemes are essentially central-difference schemes 
and have the undesirable property of being oscillatory near shocks. In 
order to stabilize them, artificial dissipation terms must be added 
which tend to smear out the discontinuities over several mesh 
intervals. Early upwind schemes include the CIR (Courant-Isaacson-Rees 
[10]) and the Godunov [11] methods. These schemes are both first order 
accurate and are not used today but were important advances and prepared 
the way for more advanced schemes. Upwind schemes are purported to be 
more physically correct since they are based on the way characteristic 
information propagates. They are also more stable near shocks and it 
has been shown that upwind difference schemes are the equivalent of 


*Numbers in brackets indicate references. 
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central difference schemes plus artificial dissipation [12]. The finite 
volume schemes popularized by Jameson [13] are also explicit central 
difference schemes which must include dissipative terms in order to 
stabilize the scheme and avoid odd-even decoupling of the mesh points. 

All of the explicit schemes are restricted to rather low CFL 
(Courant-Friedrichs-Lewy) numbers which restricts the maximum allowable 
time step they can take. This means that the number of iterations 
required to achieve a given level of convergence is larger than would be 
the case if the schemes could run at a higher CFL number. To overcome 
this limitation, implicit schemes were developed. There have been many 
such schemes developed of which one of the most significant is the Beam- 
Warming approximate factorization algorithm [14]. This scheme was a 
major advance and many later schemes were, at least in part, based on 
it. Since this scheme solves the governing equations in conservation 
form, the converged solutions satisfy the Rankine-Hugoniot jump 
relations if shocks are present. However, as is the case with all of 
the schemes discussed up to now, artificial dissipation must be present 
to stabilize the scheme near discontinuities. 

Since the numerical algorithms typically begin with an initial 
guess and then iterate towards a converged solution which must satisfy 
the specified boundary conditions, the problems to which they are 
applied are mathematically described as initial-boundary value 
problems. The discontinuities that are computed and resolved arise 
during the course of the computations and the location of these 
discontinuities is generally not known beforehand. Thus the schemes are 
frequently described as "shock capturing." Another approach that has 
proven successful is referred to a "shock fitting." In this approach, 
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the location of the discontinuity is at least approximately known 
beforehand and the appropriate mathematical relations are used to relate 
the flow conditions an each side of the discontinuity, thus providing a 
much better resolution of the discontinuity. In this study, both shock 
capturing and shock fitting methods are used and the results are 
compared. 

Separated flow has, in the past, almost always been associated with 
boundary layer separation and thus is usually regarded as a viscous 
effect. Several years ago, however, Salas [15] and others noticed that 
inviscid compressible flow past a circular cylinder computed using the 
Euler equations can separate when the free stream Mach number is greater 
than 0.4. Salas pointed out that earlier analytical investigations by 
Fraenkel [16] proved that exact solutions of the incompressible Euler 
equations, for flow past a circular cylinder, can show separation 
bubbles in front of and behind the cylinder, the size of which are 
controlled by the free stream vorticity. Salas also deduced that flow 
through a curved shock could produce sufficient vorticity to cause the 
flow to separate in some cases. He raised several questions concerning 
the validity of the computed solutions. First, are the converged 
solutions unique? Second, what is the relation between the computed 
Euler solutions to the solution of the Navier-Stokes equations, 
especially in the limit as the Reynolds number goes to Infinity? 

Kumar and Salas later compared Euler and Navier-Stokes solutions 
for supersonic shear flow past a circular cylinder [17]. The impinging 
supersonic flow contained vorticity and the separation occurred along a 
symmetry line ahead of the cylinder. The investigators found that while 
the overall size of the separation zone computed using both sets of 
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equations was similar, the internal structure was quite different which 
they attributed to viscous effects. The inviscid solution showed only 
one vortex, whereas the viscous solution showed an inner and an outer 
vortex. As the Reynolds number was increased, the inner vortex 
decreased in size and the Navier-Stokes solution became similar to the 
Euler solution. 

Barton and Pulliam also compared Euler and Navier-Stokes solutions 

for flow past a NACA 0012 airfoil at high angles of attack [18]. They 

used an implicit approximate factorization scheme which solved either 

the Euler equations or the Navier-Stokes equations in thin layer form 

with an algebraic Baldwin-Lomax turbulence model. The authors first 

describe the inviscid results for M = 0.25 and 0.4 on a coarse 249 x 

00 

41 grid and a fine 249 x 67 grid with an a = 15°. At M = 0.25 on 

oo 

the coarse grid, there was no leading edge shock but the flow separated 
and an unsteady oscillatory behavior was observed in the solution. On 
the fine grid, there was a shock but no flow separation and the solution 
converged to a steady state. At M =0.4, there was a leading edge 

oo 

shock and the solution exhibited an oscillatory separated flow which did 
not depend on the grid. Viscous and inviscid calculations were carried 
out for M = 0.301 and a = 13.5°. The Reynolds number in the viscous 
calculation was 3.91 x 10 6 . While both the inviscid and viscous 
calculations yielded a separated solution, in the inviscid case, no 
steady state solution was reached and an oscillatory behavior was noted, 
whereas the computations using the viscous equations converged to a 
steady state. The authors concluded that in this case the Euler 

solution was not a good approximation to the Navier-Stokes solution. 
They also concluded that while 171 some cases the separated flow computed 
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was the correct solution to the inviscid equations, in other cases 
numerical errors due to the use of a very coarse mesh caused the flow to 
separate and the solution was not valid. 

Separated flow on the leeside of conical delta wings at high angles 
of attack has been the subject of several investigations. Marconi 
studied supersonic conical separation using a the lambda scheme with 
shock fitting to solve the Euler equations [191. It was assumed that 
the flow is invariant in the axial direction and therefore the three- 
dimensonal equations can be solved using a two-dimensional grid. A 
shock wave emanating from the leeside of the body produced a significant 
vorticity gradient which resulted in flow separation. The separated 
flow spirals up and does not form a closed recirculation eddy as in the 
case of flow past a cylinder. Grid refinement was done without any 
significant change in the results. 

Newsome and Thomas computed and compared Euler and Navier-Stokes 
solutions for a conical delta wing at a M equal to 2.0 and a equal 
to 10° [20]. The conical assumption was also used allowing the 
calculations to be done on two-dimensional grids. The viscous solutions 
were obtained using both a central difference scheme based on 
MacCormack's explicit unsplit algorithm and also a flux vector splitting 
scheme developed by Thomas. The viscous solutions obtained with both 
schemes on a 151 x 65 grid agreed closely and showed a separated 
vortical flow on the leeside of the wing with a primary vortex due to 
leading edge separation and a secondary crossflow separation vortex. 

The Euler solutions were obtained using the same schemes after 
dropping the viscous terms and the surface no-slip boundary condition. 
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Two grids, a coarse 75 x 55 grid and the same viscous grid as before, 
were used. On the coarse grid, the MacCormack scheme solution showed a 

s 

separation vortex similar to the viscous results. Entropy was produced 
at the leading edge although there was only a small shock evident. The 
upwind flux vector splitting scheme produced solutions with no such 
separation on the coarse grid. The authors concluded that the flow 
separation obtained on the coarse grid was due to the artificial damping 
added to the central difference scheme so that the results were 
considered incorrect. Using the fine grid, both schemes produced 
similar results and the solutions showed a small vortex downstream of 
the crossflow shock. 

Chakravarthy did further calculations on the same problem using a 
TVD scheme to study the issue of inviscid separation [21]. He concluded 
that one possible explanation for the different results that Newsome and 
Thomas had obtained was the use of spatially varying time steps. He 
stated that this practice may result in unphysical transients which does 
not occur when a global time step is used. 

However, Kandil and Chuang performed a careful investigation of 
supersonic vortex dominated flows about sharp and round-edge conical 
delta wings and concluded that for the round-edge wings the damping 
coefficients used in the finite-volume code controls whether attached or 
separated solutions are obtained [22]. In addition, their computations 
indicated that the solutions obtained did not depend on whether global 
or local time stepping was used, in contrast to the hypothesis of 
Chakravarthy. 
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In view of the findings of the various investigations just 
discussed, it is clear that numerical solutions of the Euler equations 
which exhibit separated flow must be carefully evaluated to determine 
whether or not the results are realistic. This investigation therefore 
is an attempt to provide further insight into the inviscid separation 
phenomenon and specifically an attempt is made to relate the inviscid 
separation case to the viscous case. Various codes were used in the 
course of this study and the general investigation procedure is 
described in the next chapter. 


Chapter 2 


PHYSICAL MODELS AMD GRID GENERATION 

In this chapter, investigation procedures for different physical 
systems are presented and various techniques of grid generation are 
discussed. 


2.1 Investigation Procedure 

Several simple two-dimensional configurations were chosen which 
could be used to test the various schemes used to solve the Euler and 
Navier-Stokes equations for transonic flow. The first configuration is 
flow past a rearward facing step. It is known from experiment that 
incompressible flow past such a step separates at a Reynolds number, 
based on step height, less than 500 to form one or more recirculation 
vortices downstream of the step [23, 24], The incompressible flow case 
has been computed numerically by several investigators and the results 
are contained in the proceedings of a recent GAMM workshop [25]. The 
compressible case for M = 0.5 has been computed by Schmidt and 
Jameson using the Euler equations and a finite volume scheme [26]. 

The step used in this investigation differs from these cases in 
that a conformal transformation was used to generate the grid. This 
resulted in the requirement that In order to avoid singularities in the 
tranformation metrics, a sharp corner had to be avoided and instead, a 
rounded expansion corner was used. The conformal transformation was 
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used to control the location of the physical boundaries and the 
resulting suddenness of the expansion. It was found that for very 
gradual expansions the flow in most cases did not separate. Whereas, 
when the expansion corner was brought closer to a sharp corner, separa- 
tion occurred. An example of this configuration is shown in Fig. 2.1a. 

A second conformal transformation was used to produce a configura- 
tion which could be taken to represent a bump in a channel as shown in 
Fig. 2.1b. In the limit, the "bump" would become a semi-circle as will 
be discussed further in the next section on grid generation. The third 
configuration chosen for study is a circular arc inside a channel as 
shown in Fig. 2.1c. This geometry was chosen as other investigators 
have used it and it was thought desirable to be able to compare the 
results obtained in this investigation with the results obtained by 
others. Finally, a NACA 0012 airfoil with a * 0° as shown in Fig. 
2. Id was chosen as an external flow problem. 

Three computational techniques were used to solve the Euler 
equations. The first is a Gabutti shock capturing scheme which is a 
variation of the lamda scheme developed earlier by Moretti [27, 28]. 
Computer codes applied by the writer to all four test problem were 
executed on the NASA Langley computer system; some of these codes were 
vectorized to speed up the execution time. 

The second method is a shock fitting version of the first method. 
In this method, the jump conditions through the shock are explicitly 
enforced and the grid is forced to adapt to the moving shock. A 
vectorized version was developed for the rearward facing step and 
conformal bump in channel problems and used to perform grid refinement 


(a) Backward Facing Step (b) Bump in Channel 


i i { ( i < 4 SUM iU ( i i < t HU 


.. rr n rT ^ r ^ ' „„„ 

(c) Circular Arc (d) NACA 0012 Airfoil 


Fig. 2.1 Configurations Used for Study. 
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studies. The third Euler sol-ver used was an implicit flux vector 
splitting scheme which was originally developed by von Lavante with 
modifications subsequently made by the writer [47]. These schemes are 
described in subsequent chapters. 

The full Navier-Stokes equations were solved using an implicit 
upwind approximate factorization scheme developed by Rumsey [29]. The 
computer code is fully vectorized and has been shown to give accurate 
solutions for unsteady flow cases. The code was applied to test 

problems one and two as these cases showed flow separation using the 

inviscid equations. Test problems three and four showed no such 

separation and hence the viscous code was not used. 

2.2 Grid Generation 

Problems in fluid dynamics are classified as field problems since 
the solutions are represented by variables such as density, pressure, 
velocity, etc. which are functions of one or more spatial dimensions as 
well as of time. Since computers have finite memories, they can only 
solve the governing equations at a finite number of representative 

spatial locations. The equations must therefore be discretized and 
solved numerically at these locations. The purpose of grid generation 
is therefore to distribute the points at which the solution is desired 
over the spatial domain in some "optimum" sense to facilitate the 
solution. 

Usually this involves a transformation of coordinates from 
cartesian to a body fitted coordinate system. Thus the original (x,y) 
cartesian coordinates are replaced by curvilinear (s,n) coordinates' 
which wrap around the configuration to be studied. The main advantage 
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in doing this is that it is usually possible to take advantage of 
computer array data structures which is generally not the case if 
cartesian coordinates are used. However, the governing equations must 
also be transformed and additional computer memory is often required for 
the storage of metric terms and Jacobians. The transformation must be 
one-to-one and the new curvilinear coordinates are frequently assigned 
integer values such 5 = 0 , 1 , 2 , ..., I max . 

Various grid generation techniques have been devised and the ones 
which were used in this investigation include conformal mapping, 
algebraic, and elliptic methods. These will be discussed in relation to 
the specific grid generation requirements of each test problem. 

2.2.1 Conformal Mapping Technique 

For the rearward facing step test problem, the transformation from 
the computational plane to the physical plane is shown in Fig. 2.2 and 
the transformation equation is [30] 

2 * | Ck 2 - D 1/2 ♦ m (t + (c 2 - 1) 1/2 )] (2.1) 

where 

C = 5 + in 
z = x + iy . 

The metrics x„ , x , y„, y are easily obtained from Eq. (2.1) as 
5 n k n 

follows. First dz/dj; is found as 

dz/d? .1 [( t + l)/( t - 1 ) 1 1/2 

IT 


( 2 . 2 ) 
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Then x^ and y^ are found from 


(x + 1y) r = ^ * 


dz _ 1 z 
dc d? it ?-l 


Finally from the Cauchy-Riemann equations, we have 


x = y and y = - x 
5 J n n 


The Jacobian of the transformation is given by 


j . . 

a (x,y) 


5 X n y ” n x ^y 


J y_ 


5 y * - Jx , 


"k ' • J y 5 


n y Jx « 


For the conformal bump test problem, the transformation is 
Fig. 2.3 and the transformation equation is [30] 


< ■ -1 


Equation (2.4) can be solved for z in terms of z and the result 


z = i [« + fc 2 - 4) 1/2 ] 


The metric terms are then found from 


(2.3) 


(2.4) 


shown in 

(2.5) 
is 

( 2 . 6 ) 


§ = ■?[! +?/(t 2 - 4) 1/2 ] 
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For the circular arc test protrfem, the conformal transformation is a 
degenerate case of a Karman-Trefftz airfoil [31]. The transformation is 
shown in Fig. 2.4 and is 


where 


1 * 1 (c ♦ 2) + U 2 - 4) V2 

n = 2 (1 - — ) 


(2.7) 

( 2 . 8 ) 


a = 2 tan' 1 (t) 


and t is the thickness of the airfoil. The boundaries of the 

computational and physical domains are specified according to lines of 
constant £ and n. The distribution of points along the boundaries was 
specified using simple polynomial and exponential stretching 

functions. In the £ direction, a third degree polynomial was found to 
be sufficient: 

5 = C 0 + (?! - 5 0 ) Uj X + a 2 X 2 + a 3 X 3 ) (2.9) 

where £ Q and £^ are the minimum and maximum values of £ respectively, 
and a}, a 2 » and aj are coefficients which are chosen to satisfy 

X = 0 at £ * £ 

o 

X = IM-1 at £ = £i 

X = IS at £ = £ s 

d£/dx = a 1 (£ 1 -£ 0 ) at X = 0 

where IM is the number of points in the £ direction and IS is chosen 
to be X at £ = £ $ . The variable X becomes the new computational 
coordinate in the £ direction. 
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In the n coordinate direction, either a second order polynomial 
similar to Eq. (2.9) or an exponential stretching function 

" = % * (n l • "o> exp Ik ' )’ r 1 (2ao) 

was used where n 0 and nj are the minimum and maximum values of n, 
k is a stretching coefficient, and s is given by 

s = Y/(JM - 1) . (2.11) 


Y is the new computational coordinate in the n direction and JM is the 
number of points in the n direction. The polynomial stretching in the 
n direction was used to generate grids used in the inviscid calcula- 
tions and the exponential stretching was used to generate grids for the 
viscous calculations. The coefficient k was calculated using an 
iterative Newton routine to satisfy a prespecified dn/dY at n = 0 
given by 



( ni - n 0 ) k exp(ks) 
JM - 1 exp(k) - 1 


( 2 . 12 ) 


2.2.2 Algebraic Grid Generation 

Bilinear interpolation [32] was used to develop an alternate non- 
conformal grid for the circular arc test problem and the two-boundary 
technique of Smith [33] was used to develop the grid for the NACA 0012 
test problem. Both techniques require that the boundaries be initially 
specified by a distribution of points. Elliptic smoothing, as described 
in Sec. 2.2.3, was used to produce the final grids in both cases. 

The general coordinates 5 and n along the boundaries of the 
circular arc was related to the arc length along the boundaries by 
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simple polynomial and exponential stretching functions similar to those 
given in Sec. 2.2.1. The grid was made symmetric about a line passing 
normal to and through the top of the circular arc and the stretching 
function was made to satisfy the requirement that one of the lines of 
constant 5 begin at the corner of the arc and the straight lower 
boundary. 

Once the (x,y) locations of the boundary points have been 
established, bilinear interpolation can be used to locate the interior 
points. For example, the x coordinate of the interior points is given 
by 

x(? ,n ) = (1 - r) x( 0 ,n ) + r x(l,n) 

+ (1 - s) x(e , 0 ) + s x(e , 1 ) 

- (1 - r) [(1 - s) x( 0 , 0 ) + s x( 0 ,l)] 

- r [(1 - s) x(l ,0) + s x(l,l)] (2.13) 

where 5 , c, r and s all vary between 0 and 1 , and r and s are normalized 
arc lengths weighted by their relative proximities to the top and bottom 
boundaries (in the case of r) and to the left and right side boundaries 
(in the case of s). A similar equation is used to get the y coordinates 
of the interior points. Note that this method does not enforce 
orthogonality along the boundaries. However, this is later achieved 
when the elliptic technique is employed in the smoothing operation. 

The computational region for the NACA 0012 airfoil was designed for 
a "C" type grid for half the airfoil only as shown in Fig. 2.5. The 
airfoil surface is given by the following equation [34]. 
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y = 0.6 [.296? x 1/2 - .126 x - .3516 x 2 

+ .2843 x 3 - .1015 x 4 ] (2.14) 

This equation does not give y = 0 at x = 1 and it was therefore 
modified slightly and put into nested form to give y = 0 at x = 1 as 
follows 

y = .1781 x 1/2 - x [.0756 + x {.2110 - x (.1706 

- .0621 x)>] (2.15) 

The points on the boundary were distributed first from A to B usihg 
a fifth order polynomial 

r = a x 5 + a 2 c 3 + a 3 c 4 + a 4 (2.16) 

where r is the arc length and 5 the computational coordinate and aj, 

a 2 , a 3 , and a 4 are found after specifying the first derivative dr/d^ 

2 2 

at A and B, setting d r/d£ =0 at B, and requiring r to be equal to 

the total arc length from A to B at a specified value of 5 . This 

leads to a system of four equations in four unknowns which can be easily 

2 

solved. Note that since the term in 5 is not present in Eq. (2.16) 

? 2 

that d r/d? = 0 at r * 0 as well . 

An exponential stretching function was used to distribute the 
points from B to C such that dr/d? was matched at B. Along the outer 
boundary, a third order polynomial was used to distribute the points. 
Once the distribution of points along the outer boundary from D to E and 
along the airfoil surface and the symmetry line A to B to C was made, 
the two boundary technique was used to locate the points in between. 
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Referring to Fig. 2.6, the -two points a and b are first connected 
by a straight line and the points (x^, ) along it are found by linear 

interpolation 

x, - x a t t(n ) (x b - x a ) 

(2.17) 

*1 s h * tM (y b ’ y a> 

0 < t (n ) < 1 

Next, the slopes of the boundaries at a and b are computed as m b 
and and the slopes at the intermediate points are found by linear 
interpolation as m^ in the same manner as above. Straight lines passing 
through the points x^ with slopes m^ are then constructed and their 
equations are 

y s y, + nif (x - x i ) (2.18) 

Lines passing through a and b normal to the boundaries are then 
constructed which have as their equations 

y * y b -sr < x ' x b> • y s y t -5- (x - x t> (2 - 19) 

D t 

The points of intersection of these lines with the lines through x n - 
are then easily computed as x b , y b and x t , y t • Finally the grid 
points from a to b are found as 

x = x i + ci 3 fx b - x^) + o 4 (x t - x^ 

( 2 . 20 ) 

y " y i *»3 (5b ‘ y 1> + °4 (y t ‘ y 1> 

Different functions and o 4 have been experimented. One possible 


choice is 


zz 



Fi 9 . 2.6 Construction 


of Connecting Unes. 
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a 


3 



a 4 = 


2 

n 


( 2 . 21 ) 


where and are constants. 

As long as K b and K t are finite, orthogonality will be achieved at 
least on the two opposing boundaries and this technique can be extended 
to the side boundaries as well. 


2.2.3 Elliptic Smoothing 


This technique is based on the following set of Poisson's equations 


[ 6 ] 


«xx + «yy 1 P (x - y > 
n xx + "yy * Q (x ' y) 


( 2 . 22 ) 


where § and n are the computational coordinates and P and Q are the 
source terms that control the resulting grid. Since in most cases x and 
y as functions of £ and n are required, the above equations are 
transformed to 


O X - 2b X + y x = - (p X + Q X )/J 2 

55 5n no 5 n 

a y - Zb y_ + y y nn = - (P y + Q y)/J 2 

55 5n nn 5 n 

where the Jacobian J = x y - x y„ and 

5 n n 5 


(2.23) 


a 


2 


x 


n 



6 

Y 


X c * + y F y n 

5 n 5 n 


2 , 2 
x ? * y c 
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Note that If the grid is orthogonal then 

7 £ • 7 n = 0 


which implies that e is zero. 


A transformation introduced by Middlecoff and Thomas replaces P and 
Q by two new function <t> and ^ as follows [35]: 


♦ 

♦ 


t l 1 ) 2 

a 


0 # 1,2 

Y 


- ^ (-j) 


(2.24) 


with this transformation, Eqs. (2.19) become 


a x - 20 

X 

+ Y X 

= a <f> x r 

+ Y <|> X 

££ 

£n 

nn 

£ 

y n 

■ y « ’ 28 

y «n 

+ Y 

nn 

= a ♦ y 5 

+ y * y n 


The terms on the right are commonly called "source terms" and they are 
related to the mesh orthogonality and spacing. Equations (2.21) are 
solved by a SLOR (single line over-relaxation) technique which sweeps 
alternatively In each of the £ and n directions. The source terms <j> 
and i|i are initially set to zero and are then slowly changed to achieve 
the required angles and spacing of the grid along the boundaries. 

The function $ along the boundaries n = constant is varied 
according to 

4> n+1 * <j> n + Cj tane (2.26) 

where Cj is some constant and 0 is the angle at which the lines of 
constant 5 intersect the boundary as shown in Fig. 2.7. Along the 
boundaries £ = constant, the control function $ is varied according 
to the spacing ds as shown in Fig. 2.8 as given by Eq. (2.27). 


25 


4> n+1 * $ n -rc 2 (ds - ds * ) ( 2 . 27 ) 

where ds' is the desired spacing and ds is the actual spacing. 

The source term is controlled along the boundaries 5 = constant 
using a procedure similar to $ along the n s contant boundary. Along 
the boundaries n * constant, ^ could be varied to obtain the required 
spacing in a similar fashion as $ along 5 ■ constant boundaries. 


In the interior, 41 and 4/ are determined by bilinear interpolation 
similar to Eq. ( 2 . 13 ). In this case the variables r and s are related 
to 5 and n by the following third order polynomial blending functions. 


r * (3 - 25) e 2 
s - (3 - 2 n ) n 2 


( 2 . 28 ) 


This insures that 41 and ij» blend smoothly from the boundaries into 
the interior. 




Chapter 3 


GOVERNING EQUATIONS OF FLUID DYNAMICS 
3.1 Introduction 

In this chapter, the complete governing equations of fluid dynamics 
in various forms as will be used later are derived. While such a 
derivation can be found in various books on the subject, the 
presentation which is given usually depends, to some extent, on the 
author's background and orientation. For example, the classic book by 
Batchelor, which gives an excellent derivation of the Navier-Stokes 
equations, barely discusses compressible flows [36]. Furthermore, 

advanced topics such as turbulence or weak solutions to the Euler 
equations are covered in only sketchy detail in most of the introductory 
texts currently in use. This undoubtedly reflects the fact that fluid 
dynamics is an extremely broad and interdisciplinary subject which is 
both highly mathematical and also has applications in many different 
areas. 

The equations which are used in this study are the compressible 
Euler and Navier-Stokes equations in both conservative and non- 
conservative form. The viscous equations which are used are for laminar 
flow and, therefore, the topic of turbulence is not discussed herein. 
The assumptions which are made are discussed in the appropriate 
sections. The assumption that the fluid behaves as a continuem is made 
throughout. 
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3.2 Mass Conservation Equation 

The principal of mass conservation states that mass is neither 
created or destroyed. Therefore, a given differential "fluid element" 
of density p and volume 6V moving along with the fluid can be described 
from the Lagrangian point of view by the equation 

(p 6V) = 0 (3.1) 

since the mass of the fluid element remains constant. It is now more 
customary to derive the governing equations from the Eulerian point of 
view using a control volume fixed in space. Consider a control volume 
6 V with a total mass given by 

/ p dV = mass inside 5 V 
v 

The rate of change of the mass inside 6V is then related to the 
integrated mass flux through the boundaries of 6V which is 

/ ^ • n dA = net mass flux 
p 

and therefore the mass conservation law in 

It /, P <JV * - J A o * • n dA (3.2) 

Equation (3.2) is the integral form of the mass conservation law and 
applies even when p and "V are not differentiable. 

Using the divergence theorem and assuming the control volume to be 
fixed, this can be written as 

J v [jf + <H* (» *)] dV - 0 
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Since this equation must hold — for all V, it can be written in 
differential form as 

P t + ? . (p = 0 (3.3) 

It can be shown that Eq. (3.3) is equivalent to Eq. (3.1), for example 
see Karamcheti [31]. In two-dimensions, Eq. (3.3) can also be written 
as 

p t + ( pu >x + ( pV )y = 0 (3.4) 

3.3 Momentum Equation 

The principal of conservation of momentum follows from Newton's 
second law which is valid for non-relatavistic masses and states that 

t - (m h (3.5) 

where t is the applied force, m is the mass and t) is the velocity. 
For a moving infinitesimal fluid element this can be written as 

t ^ ( p * « v ) (3.6) 

where D/Dt is the substantial (or material) derivative and <sV is the 
elemental volume 

The force in the left side of Eq. (3.6) is usually considered to 
consist of body forces which are the result of gravitational or magnetic 
fields and surface forces which only act on surfaces and which are 
pressure and viscous stresses. Normally in aerodynamic analyses the 
body forces are neglected since they are negligible due to the fact that 
a moving body of air has a much larger kinetic energy than a potential 


energy. 
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It can be shown that the surface forces on a fluid elenent can be 
collectively and completely described by a second order stress tensor 
o.. (The derivation of this fact is contained in most fluid dynamic 

' J 

textbooks and will not be repeated herein). The interpretation of the 
individual components of this tensor is that 0 ^ is the force in the 
ith direction on an element of area whose normal is in the jth 
direction. It can also be shown that the stress tensor is symmetric so 
that o.. = o... This is due to the fact that an infinitesimal fluid 

' J J 1 

element cannot support moment forces as the volume goes to zero faster 
than the rotational forces which would otherwise lead to infinite 
rotational moment forces per unit volume. When referred to the 

principal axes, it is found that the off-diagonal components of the 
stress tensor are zero and that the sum of the diagonal elements, 
referred to as the principal stresses, is an invariant sun under changes 
of direction of the orthogonal axis of reference. 

The stress tensor in reference to the principal axes can be split 
and expressed as the sum of two tensors 

11 1 

22 " I T ii l 

°33 ‘ 1 T ii 

where o^, o 22 , and o 33 are the principal stresses. The first part 
of Eq. (3.7) is an isotropic tensor and the second part is referred to 
as the deviatoric stress tensor. In a fluid at rest, all of the 
components of the deviatoric stress are zero so the isotropic stress is 
- simply due to hydrostatic pressure and 
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where p is the uniform hydrostatic pressure. In a moving fluid, the 
isotropic stress tensor is still considered to be due to hydrostatic 
pressure and the components of the deviatoric sress are nonzero. Thus, 
in general, the stress tensor can be written as 


°ij * - p{ ij +t u (3 - 8) 

where is the deviatoric stress. The deviatoric stress is related to 
the motion of the fluid and in particular to the local velocity 
gradients 3U-/3X.. If it is assumed that t i • is a linear function of 

I J J 

3 u../ 3 Xj, then 

T,j = A jjkt 3^/3*, (3.9) 


where A., is a fourth order tensor coefficient. The tensor of the 
local velocity gradient can also be written as the sum of a symmetrical 
tensor, called the rate-of-strain tensor and an anti -symmetrical tensor 
which represents pure rotation. Thus 


3Uj 

3X„ 


1 3U. 3 U , 3U. 3 U 

1 r-Jl + -±) + 1 r_Ji . __ L-) 

2 '•ax a x. } ? 'a* a x. J 


3 X. 


3 X. 


e k* + 5 k* 


e k& “ 7 6 k*m “m 


(3.10) 


where w m is the angular vorticity. Thus 


ij 


A i jkjt 


(e 


1 


kfc ” 2 e ktm “iJ ’ 


(3.11) 


If is also assumed to be an istropic tensor, then the fluid is 

said to be "Newtonian" and A^^ can be expressed as the sum of the 


product of delta tensors 
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A ijk* " v 6 ik 5 jt * W 1 6 u 6 jk + w 2 6 ij 6 k* 


(3.12) 


where p, pj, and p 2 are scalar coefficients. It turns out that p is 
the molecular viscosity of the fluid which is frequently regarded simply 
as a constant but which is a function of the state of the fluid. Since 
the stress tensor is symmetric, it must be true that is also 


A ijk* must a1so be 


symmetric in i and j and that p » pj. 
symmetrical in k and i with the result that the vorticity term in Eq. 
(3.11) must be zero and therefore 


T ij = 2w e ij +w 2 e kk 6 ij (3 ‘ 13) 

where e fck = 7 • I/ 1 is the divergence of 1 1. 

The coefficient P2 is frequently called the second coefficient of 
viscosity or Lame's constant and given the symbol x. Since in a fluid 
at rest, it must be true that t .. is zero so that the mean normal stress 

* J 

is just equal to -p $ ... Stoke 's theorem assumes 

I J • 

X = - p (3.14) 


The stress tensor is finally given by 

a . . = -p6.. + 2p(e..-^7»1/6..) (3.15) 

ij ij ij 3 ij 

The rate of change of momentum inside a control volume eV is equal 
to the sum of the net momentum flux through the boundaries of 6V plus 

the change of momentum inside fiV due to the forces acting on 6V. The 

net change of momentum inside 6 V is 
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The net momentum flux in through the boundaries of 6V is 

J 3V p M . (ft 

where p H is a dyadic called the momentum flux tensor. The integrated 
force due to pressure and the viscous stresses is 

',v ( - p - " tT iJ • n) dA 
The momentum equation is then 

U / p ^ dV + / p 'H . d* « J (- p + t. .) • d* (3.16) 
dt J v J 3 v ^ ; 3 v k ij 

Equation (3.16) is the momentum equation in integral form and applies 
even when there are discontinuties in the fluid. By making use of the 
generalized divergence theorem and assumming that p, p, and ^ are 
sufficiently smooth, Eq. (3.16) becomes 

K tft (p ^ + 7 * (p ^ + P - Tjj)} dV = 0 (3.17) 

Equation (3.17) must apply to all parts of the control volume and can 
also be written simply as 

~ (p *) + V • (p H It + p) = 7 • t (3.18) 

Equation (3.18) is the Navier-Stokes equation in conservation form. If 
the viscous stresses are assumed to be zero, then Eq. (3.18) becomes 

Up ^) + 7 • (pH + p) = 0 (3.19) 

3t 

which is the Euler equation governing inviscid fluids. Equations (3.18) 
and (3.19) can be written in what is called non-conservation form by 
subtracting out the mass conservation equation which occurs in them. 


The result is 
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p ot“" 7 • °ij 


(3.20) 


which is a very compact form of the Navier- Stokes equations but which is 
still valid in general. 


3.4 Energy Conservation Equation 

The full three-dimensional mass and momentum conservation equations 
contain the five variables p, p, y, v, and w. Since there are only 
four equations thus far, another equation is needed to close the set and 
this is the energy conservation equation. In the chapters which follow, 
the energy equation will be used in different forms so a unified 
derivation is presented here. The energy equation is obtained by 
application of the first law of thermodynamics which relates the rate of 
change of the energy inside a control volume to the energy flux through 
the boundaries of the control volume and to the rate at which work is 
done on the fluid. 

The total energy of the flow field is the sum of the internal 
energy, kinetic energy, and potential energy. The potential energy is 
normally neglected in aerodynamic flows since the density of air is very 
low. Thus the rate of change of the total energy inside the control 
volume 6 V is 

■ar * ^ / v P(e + q a /2> ov 

where E t is the total energy, e is the internal energy and q 2 /2 is the 
kinetic energy. 

The energy flux in through the boundaries, due to convective flux 
of internal and kinetic energy, is given by 
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- / JV «(e tq 2 /2) M . <sk . 

The conductive heat flux is given by 

The heat flux q is related to the gradient of the temperature by 

Fourier's law which is 

q * - k 7 T 

where k is the thermal conductivity. In many derivations of the energy 
equation, a term is added to account for the heat added per unit mass 
such as would occur from an exothermic chemically reacting flow. This 
term is usually given by the following volume integral 

/ p qdV 
V 

where q is the rate of heat addition per unit mass. In this 
investigation, no such flows are considered and, therefore, this term is 
not included. 

The work done on the fluid inside sV is due to the stress tensor 
and is given by 

d * ■ ♦',,*] • 

Combining the above expressions, the integral form of the energy 
conservation equation is 

~/ v p(e + q 2 /2) dV = y [ -p (e + q 2 /2) ^ + k 7 T + (-p + x.^j-d* 
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Equation (3.21) is very general and applies even when there are 
di scontinuties in the flow. By making use of the divergence theorem, 
Eq. (3.21) can also be written in differential form as 

~ E t + 7 • [E t K k v T + ( p - t.j) K] = 0 (3.22) 

If the fluid is considered to be inviscid and non-conductive, then 
Ea. (3.22) reduces to 


it E t * 7 • [(E t tp >* 3 ■ 0 


(3.23) 


Equation (3.23) can be expanded to 

2 2 2 
3 (e + q /2) , , q , 3p , , q , . , + lS 

p - * (e . r ) — + (e + . ( p *) 


+ pV*7(e+-|-)+v-(p^) 


the second and third terms can be dropped since they Include the mass 
conservation equation to yield 


D (e + q /2) _ 
p Dt +7 


(pi/) = 0 


(3.24) 


Equation (3.24) can be split into two equations as follows. First, Eq. 
(3.20) can be written for an invicid fluid as 

dK . 

p _ + V p = 0 

This equation can be dotted with to form a scaler equation 

p V * Dt + V * vp = 0 (3.25) 

Equation (3.25) can be subtracted from Eq. (3.24) to yield 

De t. 

Dt tp7 - “ ° 


P 


(3.26) 
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which expresses the rate of change of internal energy of a moving fluid 
element. Equation (3.26) can be put in terms of the enthalpy using the 
definition of enthalpy 

h = e + pv = e + p/p (3.27) 

Substituting Eq. (3.27) in Eq. (3.26) and by making use of the mass 
conservation equation, it can be shown that the following equation 
holds. 

■>W + T& * 0 < 3 - 28 > 

Furthermore, by adding together Eqs. (3.28) and (3.24), subtracting Eq. 
(3.26), and by making use of the defintion of total enthalpy, namely, 
h 0 = h + q^/2, we obtain the following energy equation in terms of the 
total enthalpy. 

Dh„ y n 

"l!r + K - 0 (3.29) 

If the flow is steady, Eq. (3.29) reduces to 



Equation (3.30) expresses the fact that the total enthalpy along a fluid 
streamline is constant for invicid flows. 

3.5 Two-Dimensional Form of Equations 

Since the problems studied in this investigation were two- 
dimensional in nature, this section contains a summary of the governing 
equations in two dimensions and in various forms. The most compact way 
of expressing the Navier Stokes equations is in the vector conservation 
form given below 


38 


(L + F + (T = R + S„ 
t x y x y 


(3.31) 


where Q is the vector of the conserved variables given as 


Q = [p , pu, P v, e] 


and F and G are the convective flux vectors given as 


F = 


pu 

2 

p u + p 

pUV 

(e + p) u 


kT. 


G = 


pv 

pUV 

pv 2 + p 
(e + p) v - kT 


and R and S are the viscous flux terns which are 


R = 


xx 


xy 

Ut + vt 
xx xy 


S = 


xy 


yy 


^xy + ^yy 


The viscous stress terms in Eq. (3.31) are 


xx 


2 f , 3 u 3 V-v 

( Z Fx ' ay) 
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yy 


2 — , 0 3 v 3 iu 

7 U ( 2 iy “ ax-J 


(3.32) 


xy 


- u 


rHUi* 

''Sy ax - 1 


The total energy per unit mass is the sun of the internal energy and the 
kinetic energy as follows: 


e 


= i + 


2 2 
u + v 


(3.33) 


Equations (3.31), (3.32) and (3.33) together contain five scalar 
equations in the nine quantities p, u, v, p, e, i, t, k, and p. In the 
model being used, the thermal conductivity k and the molecular viscosity 
u are considered constant so two additional equations are needed to 
close the system. The first is the perfect gas equation of state 


P * p R T 


(3.34) 


where R is the gas constant. If the fluid is assumed to be calorically 
perfect, then 


1 = c v T 


(3.35) 


and consquently by Eq. (3.34) and the relation 


v y - 


(3.36) 


it holds that 


i = p/[(y-1)pL 


(3.37) 


The governing equations can b£ nondimensionalized by referring them 
to suitable reference quantities as follows: 


= x/*.. 


y’ = y/i, 


t * = t 


a o /l o 


1 
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r 



u' 

= u/> 0 


V' -v/a 0 p' = p/p 0 

(3.38) 



p' 

■ P/fco a o> 


e ' s e/( P 0 a o J u ' = v/v o 




r 

- I/'I 


V - T/T 0 


c 

' 

The resulting equations 

after nondimensionalizing are: 



r ~ 


ia; + 

FT 

3 F 1 
3 X 1 

+ 3 G* _ l r a R ’ as\ 

37" 'Re lax" 37 rJ 

(3.39) 



where 

Q’ 


Lp , p u , p v , e J 
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and 



2 


1 (2 


3 u' 

ifx T 


3_V 

3y 


•) 


T 


yy 


2 , 


(2 


3 V 1 

W 


a_u 

3 X 


■) 


x 1 = r iX + K-) 

yy '-Fy 3x • 1 

The reference Reynolds number is gi.ven by 


Re. 


a o Po l o 


The equation for the total energy, Eq. (3.33), becomes, after nondimen- 
sional izing, 


. 1 2 . „i 2 


e- , i* 


(3.40) 


and the equation of state, Eq. (3.37), becomes 


i 1 = p 7 C (y - 1 ) p’] 


(3.41) 


Thus the original equations are replaced by equations which have the 
same form but which now include the reference Reynolds number in the 
momentum equations and the reference Peclet number (Re 0 Pr 0 ) in the 
energy equation. 

The Euler equations for inviscid fluids are obtained from Eqs. 
(3.31) and (3.39) by dropping the viscous flux terms R and S. Thus we 
have 


Q t + F x + G y 


0 


(3.42) 
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Equation (3.42) is written out as~ 

p t + ( pu) x + ^ pV ^y = 0 

(pu) t + (pu 2 + p) x + ( P uv) y = 0 

( P v) t + ( P uv) x + ( p v 2 + p) y * 0 (3.43) 

e t + [(e + p)u - kT] x + [(e + p)v - kT] y = 0 

Eauation (3.43) are in conservation form and include an energy equation 
in terms of the internal energy. In the three schemes used to solve the 
Euler equations, this form of the energy equation was not used. 

Instead, energy equations in terms of either the enthalpy or entropy 
were used. 

In the following discussion and derivation, the primes are dropped 
for convenience and it should be understood that all quantities are non- 
dimensional . 

The next step is to transform the equations to general curvilinear 
coordinates in conservation law form. Equation (3.39) is transformed to 

Q+F+G = R + (3.44) 

■t 5 n £ n 

where 

Q - Q/J 

F - (F £ x + G £ y )/J 

*6 = ( F n + G n )/J 
x y 

R - (R? x ♦ Sc y )/J 
*5 = (Rn x + Sn y )/J 
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Chapter 4 

SCHEME 1 - SHOCK CAPTURING GABUTTI 
4.1 Introduction 

The first computational scheme used in this study is the Gabutti 
scheme which is a refinement of the x-scheme developed earlier by 
Moretti [28], The scheme solves the time dependent compressible Euler 
equations in non conservation form by an explicit predictor-corrector 
method 1 . Gabutti 's method improves upon the x-scheme by extending the 
stability range substantially and also by enabling it to satisfy the 
shift condition for a CFL of one which is 

n+1 n 

u i = u i-l * 

The scheme has good shock capturing properties in that the discontinuity 
is usually spread over no more than two or three mesh points. The jumps 
in density, pressure, and velocity through the shock are not correct due 
to the nonconservative nature of the scheme, however, and shock fitting 
must be used to get the correct jump relations. The scheme solves the 
mass conservation, and x and y momentum equations in non-conservative 
form as 


p t + Up x + vp y + p < u x + v y) = 0 

u t * uu x * vu y ■ ° (4 - 1) 
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*t ♦ uv x * vv y -> i. P y - 0 . 

If the flow is assumed isentropic along streamlines, and if we use the 
change of variables P = in p, then the above equations can be written as 


P 

u 

v 


t 

t 

t 


+ 

up 

+ 

vP 

+ Y (u + 


X 


y 

T X 





2 

+ 

uu 

+ 

vu 

+ p 


X 


y 

Y X 





,2 

+ 

uv 

+ 

vv„ 

+ — p 


X 


y 

y y 


0 

0 

0 . 


(4.2) 


The energy equation in its usual differential form is not used, but 
rather the assumption is made that the flow is isenthalpic. In this 
case, the square of the speed of sound a^ is related to the velocities 
by the steady state energy equation 



+ 



(4.3) 


where h Q is the total enthalpy. If the pressure p and the density p are 
nondimensional ized by dividing by their stagnation values p 0 and p Q , 
then 

a* 

h o * ft ' T^r • < 4 ‘ 4 > 


4.2 Transformation of Equations 

In order to use the governing equations, it is first desirable to 
transform them to (c,n) coordinates. Using Eqs. (2.4) and the chain 
rule, the mass conservation equation becomes 

P. + UP + VP + Y [U + V 1 = P c (4.6)' 

t 5 n 5 n s 

where U and V are the contravariant velocities 
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U 

■ f~x u + V 

(4.7a) 


V 

* V + v 

(4.7b) 

and 

P s 

* -1 J tg 4 u + g 5 v] . 

(4.7c) 


The terms g 4 and g 5 are two of 11 transformation terms that appear in 
the governing equations after transformation and these are given by 


g l 

= 

Z . 2 

g 2 

= 

n n 

g 3 

= 

X ti + y c 

^4 

= 

X 5E • X Cn y « ’ y 

g 5 

= 

X £ n y 5 - X nn y E ' y 

g 6 

= 

y n X S 5 ‘ x n y « 

g 7 

= 

y x - x y 
n nn n nn 

g 8 

= 

2 (y x - x y ) 
n 5n n £n 

g 9 

= 

x y - y x 

10 

s 

X 5 y nn - y 5 

11 

* 

2 (X E y 5n - y 5 V 


(4.8) 

The transformation of the two momentum equations is somewhat more 
complicated. First, they are written in vector form as 


Q t + uQ x + vQ y + S = 0 


(4.9) 
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where 

The vector Q can also be written as 


Q - [u, v] T and-S * ^ CP X . P y ] T . 


where 


Q = 


» 

u 


X X 


‘U 


- 

5 n 



V 


y e y 


V 

* ■ 


L 5 n J 


m 


* T Q 


T = 


Q « 


x x 

5 n 

y_ y 

5 n 

U 


(4.10) 


a a ~-l 

Multiplying Eq. (4.10) by T gives 

PW <V« 1 ,ow PV 4 »v /V 1 

Q t + T [u (T Q) x + v (T Q) ] + T S = 0 . 

Note it is assumed that T is independent of time. The previous 

equation is equivalent to 

Q* + T [U(T Q) + V(T Q) ] + T S = 0 

t 5 n 

which can be manipulated to get 

Q+UQ +VQ + $ = 0 (4.11) 

t 5 n 

where 

PS* 1 _ M IV • I 

S = T [U T OUT 0] + T S . 

5 n 


Equation (4.11) written out is 

u t . uu 5 + vu n ♦ ^ [g 2 P 5 - g 3 p„3 - u s 

2 

v t . m % ♦ w ♦^-[ 9 1 P n - 9 j P 5 ] - V s 


(4.12a) 


(4.12b) 
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where 

U $ * -J (g 6 U 2 + g, V 2 + g g UV) 

(4.12c) 


v s » -J (g 9 u 2 + g 1Q V 2 + g n UV) . 

(4.1 2d ) 


If the transformation from (x,y) to (c,n) is conformal, then the above 
two equations can be reduced using the Cauchy-Riemann equations to 

2 2 2 

IL + UU + VU + — P - . u — y- L J. - UV J, = 0 (4.13a) 

t £ n y 5 Z1 Z 


V. + UV + VV + — P - UV J 
t {■ n y n 




0 


(4.13b) 


where 

Jj = J^/J 

J 9 = J /J . 
2 n 


The x-scheme introduced by Moretti is next used to put Eqs. (4.6) 
and (4.12) into a quasi-characteristic form. First, the time derivative 
of P is split into two parts and P’Ji such that 

f\ * u P { * y u { - \ % (4.14) 

■ 7% • <«•«> 

Equations (4.12a) and (4.14) can be written as 


where 


^ 5 


A 



y 

u 



[ (Ja) z 




(4.16) 
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and S} is the vector of the remaTning terms put on the right side of Eq. 
(4.16). The eigenvalues of matrix A are 

± , , ** 

X = U t . 

Equation (4.16) can be diagonalized as 

+ [S] Cx] [S]" 1 Q* = S. (4.17) 

t 51 

where [S] is the matrix of the eigenvectors. If [x] is split into its 
positive and negative eigenvalues, Eq. (4.17) becomes 

35 ♦ [S] [X + : [S]' 1 + [s] Cx'] [s]' 1 Q 5 = S . (4.18) 

t s ^ 1 

Similary, Eqs. (4.12b) and (4.15) can be written as 


where 


v t *ttr - s 

t n 2 




■ [(J ») 2 




(4.19) 


The eigenvalues of Eq. (4.19) are 

Q i = V ± a 2 . 

Equation (4.19) can be diagonalized and split into positive and negative 
components. The result is 

Q? * [R] Co + ] [R ]' 1 V * [R] On'] CR ]' 1 ? ■ S, . 

t n n 2 


(4.20) 



Equations (4.18) and (4.20) can t>e added together to get the final set 
of equations which are given below 
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* + i [xV + x'p; + nV + n’P“ + rV (x + u 

t 2 u E E n n aJ« I 


5 

+..+ 


V 
c 


+ Or V „ - «'0 

au^ n n 


= P. 


x"u-)3 


u t + \ [x\ + + x'u; + (X + P p + - x‘p;)] 


5 5 Y 

+ VU - -M- g~ P 
n y 3 n 


£ 

U. 


V. + i [nV + q‘v” + (n + P + - n‘P')] 

t 2 *- n n y n n 


+ UV - -M- g. P 
n y 3 5 


= V. 


(4.21) 


where P $ , U s and V s are unchanged but 




n n 


4.3 Discretization of Equations 

Equation (4.22) is now in a form that allows a proper discretiza- 
tion of the spatial derivatives according to the signs of the eigen- 
values multiplying them by taking into account the domain of dependence. 
In two-dimensional unsteady flow, the governing Euler equations are 
hyperbolic in time and the solution at any point X. at time t + At 
depends on the solution at time t within the area dA as shown in Fig. 
4.1. The conoid from the solution surface at time t to point X at 
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t + At is called a Monge cone— and it depicts the way the solution 
propagates in time [37]. The numerical algorithm should model the way 
information propagates which means it must take into account the proper 
domain of dependence. Therefore, the information used to update each 
mesh point should be taken from adjacent mesh points in a manner that 
reflects the physical propagation of information. In supersonic flow, 
for example, pressure waves do not propagate upstream and no information 
from downstream points should be used to update the solution at a given 
mesh point. 

The Courant-Frederick-Lewy, or CFL, stability condition reflects 
the requirement that the region of dependence must be at least as large 
as the analytical domain of dependence [6]. In addition, Moretti also 
refers to the "law of forbidden signals" and states that in addition to 
satisfying the CFL stability condition, a numerical scheme should also 
satisfy this law [28]. llhat this means in practice is that in 
supersonic flow, the information that is used to generate the updated 
solution at a particular mesh point should not come from points 
downstream of it. This is not always possible but the present scheme 
does attempt to follow this principal as closely as possible. 

Equations (4.21) are solved using the three-step predictor- 
corrector scheme [27]. The scheme will be described using the one- 
dimensional linear wave equation as a model problem 

u t + a u x = 0 . (4.22) 

In step 1 of the scheme, the spatial derivative u x is evaluated by 
either a two point backward or forward difference, depending on the sign 
of the characteristics as 


5 ] 



Fig. 4.1 Characteristic Monge Cone. 
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U i ' u i-l 

u x = 7u i = —IT- if a < 0 

u i+i " u i 

u « AU. = if a > 0 . (4.23) 

A ! Zi A 

A predicted value is then calculated using the u t that results from 
Eq. (4.22). For example, if the wave speed a is positive then 

u. = uj + u t At (4.24) 

where u t = “ a u x 

In step two of the scheme, the values of u^ at time level n are 
used again to compute u t . Depending on the wave speed a, u x is computed 
as 


2u 1 - Kl + u i-2 

AX 


$ 3 > 0 « 


(4.25) 


+ C2 

AX 


, a < 0. 


(4.26) 


Then, compute u t as 


u t = - a u x 

In step three of the scheme, the predicted values u.j are used to 
compute a predicted u^ with two step backward or forward differences, 
depending as in Eq. (4.23) again on the sign of the characteristic speed 
a. The final update is then made using both the u t computed in step two 
and u t as follows 

u" +1 = uj + ! (u t + u t ) * (4.27 ) 
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In this example, the sign of the- wave speed a determines the direction 
in which the spatial derivatives are taken. In Eqs. (4.21), the 
characteristics x + and \~ , and n + and a" perform this function. The 
scheme satisfies the so-called "shift" condition for v = 1 which is 

u n+1 = u n 
u i u i-l ‘ 

Thus, for v * 1, the scheme properly convects waves along the 
characteristic dx/dt = a. 


4.4 Stability Analyses 

The stability of the scheme was analyzed using the classical von 
Neumann method [38]. In this method, the error is expressed as a 
Fourier series and the growth of the error in tine is examined. The 
total error in a numerical calculation consists of both discretization 
error and round-off error. The discretization error arises from the 
fact that what one is actually solving in a numerical calculation is not 
a differential equation but a difference equation. If U(x,t) is taken 
to be the exact solution to the governing PDE and if u(i'ax, nat) is 
the solution to the approximating difference equation carried out to 
infinite precision, then (U - u) is the discretization error. If this 
difference goes to zero as ax goes to zero, then the difference 
equation is said to be consistent with the PDE. 

Roundoff error arises from the fact that the calculations cannot, 
in practice, be carried out with infinite precision and the calculations 
must be rounded-off to some finite number of decimal places. If N is 
called the actual finite precision numerical solution, then (u - N) is 
the roundoff error and U - N is the sum of the discretization error and 
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the round-off error. The numerical calculations will be stable if the 
inevitable errors which are introduced are damped out (i.e. decay with 
time) . 

Although the von Neumann stability analyses is not the most 
rigorous method, since it ignores the boundary conditions, it is used 
frequently because one is usually interested in the stability of the 
basic scheme as boundary conditions may change. To examine the effect 
of the boundary conditions, the matrix method may be used [39], In this 
method, it is necessary to define an amplification matrix A as 

u n+1 - A u u 

where u n+1 and u n are the vectors of the solution. Then the 
eigenvalues of A are examined, and for stability, it is necessary that 
the modulus of all the eigenvalues be less than one. 

Stability analyses of the scheme were done in both one- and two- 
dimensions. For the one-dimensional analyses, the model Eq. (4.22) was 
used. For the two-dimensional stability analyses, the following model 
equation was used. 


u t + a u + b u = 0 (4.28) 

It turns out that the maximum CFL number for which the scheme is stable 
is two for the one-dimensional equation and one for the two-dimensional 
model equation. The steps in the analyses will only be shown for the 
two-dimensional analyses. 
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The first step in the analyses is to combine the three steps of the 

Gabutti scheme into a single step so that u? + J Is expressed in terms 

' *J 

of the u.j j at time step n. When this is done, the result is 


u 


n+1 

i,J 


n 1 r / « n 
u ij -7M 3 u ij 



+ 



*8(3 u" j 






2 a - u Mjj - * u 1 _ 1< j. 1 J] (4-29) 


where 


a = a AX/&t 


6 = b Ay/At 

If e" is assumed to be an initial error distribution and N? is the 
solution to the difference equation which satisfies Eq. (4.29) exactly, 
then the error at time level n+1 must also satisfy Eq. (2.29). Thus, 
the error will not grow provided the solution to the difference equation 
is stable and bounded. 

Next, assume that 

u" = r n exp [I (k x i ax + k y j Ay)] (4.30) 

I = vM 

is a periodic representation of the solution with r n as the amplitude, 

k y and k„ as the wavenumbers, and 9„ = k„ a„ and e w = k Ay as the 

A j x x x y y 

phase angles. Substituting Eq. (4.30) into Eq. (4.29), is found that 
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r * |l - \ [a (3 - 3~exp(-Ie x ) + exp(-2Ie x ) ) 

+ 6(3 - 3 exp(-Ie ) + exp(-2Ie )) 

+ a 2 (l - 2 exp(-Ie x ) + exp(2Ie x )) 

+ 0 2 (1 - 2 exp(-Ie ) + exp(2Ie y )) 

+ a 6(2 - 2 exp(-l 0 x ) - 2 exp(-Ie y ) 

+ exp ( - I (© x + 0 y )))]| . (4.31) 

Define the amplification factor as 

uj +1 - G uj (4.32) 

from which it Is apparent that G is the same as r. The maximum 
amplification factor G for values of o and 0 ranging from .5 to 1.5 
was found by solving Eq. (4.31) on a computer for values of o and 0 
ranging from 0 to 360 degrees. The results are shown in Table 4.1. 

From this, it is seen that the scheme has a maximum CFL number for 
stability of one in two-dimensions. 



Table 4.1 Results from Stability Analyses 


a 

B 

G 

.5 

.5 

1.0 

.5 

1.0 

1.0 

.5 

1.5 

1.0 

1.0 

1.0 

1.0 

1.0 

1.5 

3.5 


1.5 


1.5 


7.0 



Chapter 5 


SCHEME 2 - SHOCK FITTING GABUTTI 
5.1 Introduction 

Shock fitting, in contrast to shock capturing, does not attempt to 
apply finite differencing across shock waves and instead imposes the 
correct Rankine-Hugoniot jump relations at the discontinuities, the 
location of which is assumed known. Shock fitting methods seem to have 
evolved alongside and in conjunction with shock capturing methods. 
Shock fitting is still frequently used in supersonic problems where a 
bow shock wave develops ahead of the body. The flow ahead of the shock 
is usually taken to be freestream uniform flow and so there is no need 
to apply the time-dependent method to this region. The jump relations 
are applied at the shock, which is fitted as a computational boundary. 
If all of the shocks in the flow field are fitted, then it can be argued 
that the governing equations can be solved numerically in either 
conservation or nonconservation form in the smooth regions of the flow 
field. 

Thus, the use of shock fitting can lead to a reduction in computer 

time from (1) the elimination of the need to apply finite differencing 

to the freestream uniform flow ahead of the bow shock in supersonic 

problems and (2) the use of the primitive variable form of the governing 

equations which may be somewhat less expensive to solve numerically. 

However, the subsequent development of implicit methods has, to some 
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degree, removed the advantage of-solvina the non-conservative equations 
because ultimately one Is after the solution to the problem In the least 
amount of computer time and convergence rates are perhaps more relevant 
to this quest than the cost to do a single iteration. 

In this study, the shock which develops, is embedded in the flow 
field and shock fitting is used on a grid which has the shock aligned 
with one of the lines of constant X (in the computational domain). The 
shock is allowed to move and adjust its position and characteristics to 
the evolving flow field. Recently, Moretti has shown that it is 
possible to do shock fitting on a grid which is stationary such that the 
shock is not aligned with any of the mesh lines [40]. This method has 
also been referred to as "front tracking" and is used to describe 
methods for handling discontinuities such as weather fronts and oil flow 
in porous media [41]. 

The shock fitting scheme developed for use in this investigation 
treats the embedded shock which develops as a discontinuity aligned with 
a line of constant X and the Rankine-Hugoniot relations are used to 
relate the upstream and downstream flows through the shock. 

This scheme is identical to the scheme described in Chap. 4 with 
five exceptions. First, the initial conditions are the converged 
results of a previous run of scheme 1. Second, the grid is aligned with 
the embedded shock that forms in the initial run. Third, the governing 
equations along with the Rankine-Hugoniot relations are used to 
calculate the shock velocity and acceleration at each point on the shock 
front and the grid is dynamically adapted as the shock changes its 
position. Fourth, the Rankin-Hugoniot relations are used to calculate 
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the jumps in pressure, density, velocity and entropy through the shock, 
and fifth, the equation describing the convection of entropy is included 
with the set of equations to be solved 

I ■ vv ,s , ■ 0 • < 5 -" 

When Eq. (5.1) is transformed to (s,n) coordinates, it becomes 

S. + US + VS =0 (5.2) 

t ? n 

where U and V are the contravarient velocity components. The entropy is 
related to the nondimensional pressure and density by Eq. (5.3) 

S = tn p • y in p . (5.3) 

The unknowns in the final set of equations are the log of the pressure 
P, the velocity components U and V, and the entropy S. Since there are 
four equations, the system is closed. 

5.2 Transformation to Shock Fitted Coordinates 

The solution to be used as the initial condition to the shock 
fitting calculations is associated with some curvilinear coordinates £ 
and n which, in general, are not aligned with the shock. The solution 
must, therefore, be interpolated onto a new coordinate system which is 
aligned with the shock. These coordinates will be called X an Y herein 
and the transformation from the previous to the new coordinate system is 
given as 


(t, 5, n) ♦ (T, X, Y) . 



61 


The 
Y = 


transformation relations to-be used will have X = X 
Y (t,n) and, therefore, by the chain rule. 


3_ 

3 1 

3_ 

3? 

l_ 

3 n 

The system of governing 
P T * UP X * VP Y ♦ t [U x x ? 


+ x + Y 

3 T 3 X t 3 Y t 


= It x n + It v t, • 

equations which results is 


V ¥ X 
X n 


Y V Y ] 
Y n 


P 

s 


(Ja)‘ 


U T + UU V + VU v + 

T X Y y 


Cg 2 V € * g 3 ( V n + V n )] 


V T + UV X + vv y + 


(Ja) 


[g (P X + P V Y ) - g P ¥ X ] 
y 1 X n “ n 3X5 


(t, e , n ) and 


(5.4) 


= U (5.5) 
s 

= V 

s 


S T + US + vs =0 
T x y 

where 

0 = X + UX + VX 
ten 

V = Y, + VY 
t n 

and P $ , U s and V s are as previously defined in Eqs. (4.7c) and (4.12c 
and d). 

The desired coordinate transformation in the 5 direction is 

X = X ( C$ ) 

where e $ is the location of the shock which is function of n and t. 
A second degree polynomial given as 
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5 - 5 „ * <*1 * « 0 > (a l X * a 2 (5 - 6) 

was used where g and g, represent the minimum and maximum values of 
g and a^ and ^ are coefficients determined such that 


5 = 

g $ when X = X $ 

6 = 

gj when X = . 

It is easily found that 


a l 

= bj + b 2 f 

a 2 

= b 3+ b 4 f 

where 

b i 

• - X s /(XjD) 

b 2 

= X 1 /(X S D) 

b 3 

= -bj/X, 

b 4 

iH 

X 

X, 

CVJ 

-O 

1 

II 

D 

■ *1 - X s 

f 

«s -«o 
«1 "«0 ' 


Using the above expressions, the transformation Eqs. (5.8) are 
determined as 

\ - l/[(5, - «„) <aj + 2 a 2 X)] 
x * ; - (b 2 X + b 4 X 2 ) x c f n 


(5.8) 
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X t - - (b 2 “X + b 4 X 2 ) X^ f t . 

A transformation from n to Y was done for problems 1, 2, and 3. 
Since £ and n are conformal coordinates and are independent of the 
stretchings used to produce the original grid, then it is necessary to 
use some sort of transformation from n to Y if it is desired to have 
a Y = 1 . A second order stretching was used as follows 

n • n Q + (n x - n 0 ) (a 4 Y + a 5 Y 2 ) (5.9) 

where n 0 and r\^ are the minimum and maximum values of n and a^ and 
a 5 are coefficients chosen to achieve the desired stretching. In 
problem 1 , n 0 » 84 and were made to be functions of time so that the 

position of the lower boundary as well as the degree of clustering near 
the lower boundary could be changed dynamically. 

5.3 Calculation of Shock Acceleration* 

To calculate the shock acceleration, the governing equations are 
used along with the Rankine-Hugoniot relations in a fairly straight 
forward manner. First, the governing equations are written in vector 

form as 

0 t + A 0 X = R (5.10) 

0 * [P, U, V] T 


*The writer is indebted to Mr. M. D. Salas for showing how to calculate 
the shock acceleration. 



where 


U 


(Ja) / _ u _ y X 

(g, x - g X ) 

Y 2 c 3 n 

(9, X - 9, X ) 

Y In 3 5 


R * [Rj, R 2 , R 3 ] T 


YX f Y X 

C n 

0 0 


Ri ■ Pc - V P v -Y U. 


U s ■ V U Y + g 3 P Y 


(Ja) 


V s . V Vy - 8l Py . 


The matrix A has the eigenvalues 
Xj = 0 - a 
x. - 0 + a 


x „ = U 


a = 


Ja [ 8l x 2 - 2 g 3 x 5 x n ♦ g 2 X 2 ] 1 / 2 


= Ja [SQ] 1/2 


and can be diagonalized as 


A = SA S' 


where 


h sq 

9 3 \ 

9j X - 9l X 


9 Z X, 


fc Sq 0 

=2 x 5 - s 3 x „ x n 

«1 X „ ' *3 X { - X { 


s 



65 


and SQ Is the quantity in square -wots. Equation (5.10) can be put into 
characteristic form as 

S" 1 0 t + A S" 1 Q x = S” 1 R . (5.11) 

The equation corresponding to Xj is then 


P t + P X + B 2 ^ U T + x U X } + b 3 ^ V T + x V 

= + 6 2 ^2 + B 3 = P 4 (5.12) 


where 

6 2 = -TX 3 /a 

So = - y X /a 
3 n 

and P~, U‘ and are forward differences. 

Figure 5.1 shows the shock in the (e,n) coordinate system with 
unit normal and tangent vectors. The covariant base vectors are 



n n 


The contravariant base vectors are 



i + Y 
n 


j 

j . 


? - «. 




e n 


n i + n j • 
x y 


(5.13) 


(5.14) 


The unit normal at the shock is given by 



67 


M 

f ? + f j 

„ x y 

/I 2 

f + f 
x y 

f I 5 + f * n 

« -i — - — a — (5.15) 

/vf • vf 

where f = £ -5 s (n, t) = 0 (5.16) 

and the length of grad f is 

|vf I - [f 2 l 5 . e 5 + 2 f f I 5 . I n + f 2 e n . e n ] 1/2 
* h n ‘ 

The unit tangent vector is defined by 

N . f = 0 


and therefore it is found that 


T = 


- f e + f e 
n 5 n ti 


r ..2 + + _ . . ♦ + ,2 + ♦ , 
[fe.e-2ffe«e+fe»el 
1 n £ 5 5 n ? n £ n n J 


1/2 


5 s + \ 

s n * n 


(5.16) 


where r is the shock slope and h t is the square root term above. 

j * # 

n 
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The velocity relations at— the shock are obtained using these 
equations. The flowfield velocity is 

1/ = u i + v j 


= U e + V £ . (5.17) 

5 n 

Letting the contravariant velocity of the shock be U s , the velocity 
relative to the shock is 


^ = (U - U ) e + 1 / e . (5.18) 

s s ^ n 

Note that U and V are different on each side of the shock. The velocity 
component normal to the shock is 

(u-u s ) - ?s V 

V . N - 2_ . 

s h 

n 

= V n (5.19) 


The component tangent to the shock is 



5 (U - U ) + V 
*s s 


n 



(5.20) 


Since the tangential components of velocity must be equal on both sides 
of the shock, we have 


Ms +V 1 " M: 


+ V, 


(5.21) 
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where the subscripts 1 and 2 denote the upstream and downstream sides of 
the shock respectively. Equations (5.19) and (5.20) can be solved for 1) 
and V to get 


“ s * 


V, h. 

n n * s t t 
n 


(5.22) 


V = 


V t h* - r V h n 
t t s s n n 




(5.23) 


The jump in pressure through the shock is given by 

= 1 + r£r (M i- « = Pr 


f2 

Pi 


(5.24) 


where is the Mach number relative to the shock and is given by 

V . N (U 1 - u s> - 5 S V 1 


1 


a h 
1 n 


Equation (5.24) can be expressed in terms of the log of p as 

P 2 = p i + * n (P r ) • 

The derivative of P 2 with respect to time is then 

% ■ \ * r; M i \ • 


(5.25) 


(5.26) 


The tern is equal to 


C 1 + c 2 U s. 


(5.27) 
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where U $ is the shock acceleration and c^ and c 2 are coefficients 
found to be 


C 1 


U U ‘ 5 s V 1 ' *s 'I 


V, 

n t 


a l h n 






Substituting Eq. (5.27) into Eq. (5.26) and rearranging gives 


bj ♦ b 2 u 


(5.28) 


where 


* \ + T r "i c i 


’2 * 'r^r 1 H i c 2 


The jump in the normal component of velocity relative to the shock is 
given by 


^ . N 

Jl 

^ . N 


= u 


U„ - U - 5 V 0 
2 s S s 2 

n 

U. - U - c V. 
1 s *s 1 


(5.29) 


Using Eqs. (5.21) and (5.29), the contravariant component U 2 is found to 
be 


U, 


u $ * 5 S * Uj ♦ u r (U, -» s .V lU ) 

n n n 


(5.30) 
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where the jump in velocity through the shock u p is given by 

(y - 1) M? + 2 

u = . 

r (y + 1) MJ 


(5.31) 


The time derivative of Eq. (5.30) is then found to be 


where 


\ * b 3 +b 4% 


(5.32) 


= t V l «s tV l 5 s * U 1 «S t 2 U l«S «! 


t 


ri n. 


♦ c, (U, 


v s 


) ♦ 


u r' U l. 


- V, - V, r ) 




5 S ]/(l * ) 

n t n 


t >4 ■ [1 ♦ c 4 (U t - U s - Vj E s ) - u p ]/(l t 5 2 ) 

n n 

c 3 - - 4 Cj/fr + 1) M 3 ) 
c 4 = c 3 c 2 /c l * 

In a similar manner, taking the time derivative of Eq. (5.21) leads to 

the following expression for V 9 . 

C X 

V 2 t = b 5 +b 6 U s t < 5 - 33 > 


b 5 = \ * ' U 1 - V «s * (U l f • b 3> 

t n t t n 


- b 4«s • 


where 
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Finally, substituting all of these expressions into Eq. (5.12) and 

solving for U gives 
s t 


*4 ~ b l “ X P x ~ 6 2^ b 3 + x V ~ 6 3 (b 5 * x V 
b 2 + e 2 b 4 + B 3 b 6 


(5.34) 


Thus an expression has been obtained for the shock acceleration. The 
shock position is then updated using the shock acceleration to update 
the shock velocity as described next. 

5.4 Updating the Shock Position 

The method used to update the shock position was arrived at after 
extensive experimentation. It was discovered early that the calculation 
of shock acceleration and the jumps in fluid density, velocity, 
pressure, and energy through the shock are all very sensitive to the 
shock velocity and the slope of the shock with respect to the oncoming 
flow. If each point along the shock is computed independently of the 
others, wiggles and oscillations tended to develop in the shape of the 
shock which quickly destroyed the calculations. In order to prevent 
this from occurring, polynomial least squares smoothing was introduced. 

A related problem which had to be addressed was how to update the 
grid above the shock. Since no acceleration is computed for the points 
above the shock along the line of constant 5 which aligns with the 
shqck, some artificial means had to be introduced to move these points 
so the grid would not have a discontinuity in it. 

The first step in updating the shock position is to smooth the 
computed shock accelerations computed from J = 1 to J = JS, the last 
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shock point. This was done by a -weighted least squares approximation to 
the function values at the JS points using orthogonal polynomials [42]. 
Then, using the smoothed values of the shock acceleration, the updated 
shock velocities and positions were computed in a three-step predictor- 
corrector manner analogous to the Gabutti scheme. 

In the first step, the following loop is executed 


Do from 

J = 1 to 

J = JS 



XSTN(J) 

= XST(J) + 

XSTT(J) * At 


XSN(J) 

= XS(J) + 

XSTN(J) * At 


where XS, XST, and XSTT are the shock position, velocity, and the 
acceleration at the beginning of the time step and XSN and XSTN are the 
predicted values of the shock position and acceleration. The predicted 
values of XSN are then smoothed using orthogonal polynomials up to 
degree three as shown in Fig. 5.2a. 

The points above the shock along the connecting line from J = JS + 
1 to J = JMAX must now be computed and this was done as follows. First, 
XSN at the top boundary (J = JM) was determine as 

XSN(JM) = .9 * XSO(JM) + .1 * XSN(JS) . (5.36) 

This insures that eventually XS(JM) will equal XS(JS) and that the top 
end of the connecting line will not move too rapidly. Next, the points 
in between XSN(JS) and XSN(JH) are calculated as 

XSN ( J ) = (1 - 3H 2 -H 3 ) * XSN(JS) 


+ (3H ? - 2H 3 ) + XSN(JM) 


(5.37) 
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where ~ 

H ' ^ n J " n J$^ n JM ‘ n JS } * 

This insures that the slope of the connecting curve £ s is zero at J = 

n 

JS and J = JM as shown in Fig. 5.2b. Finally, the points XSN from J = 1 
to J = JS are smoothed to blend the connecting curve with the shock as 
shown in Fig. 5.2c. The smoothing this tine is done by orthogonal 
polynomials up to degree four. 

In the second step of the Gabutti scheme, the shock accelerations 
are again calculated at each point then smoothed by second degree 
orthogonal polynomials and stored as XSTTS. The predicted values of the 
shock position and the points along the connecting curve are now 
substituted for the original values of XS. 

In the third step, the same procedure is followed as in the first 
step except the calculated values of the shock acceleration are first 
smoothed and then combined with the predicted values from step two as 
follows 

XSTT(J) = j (XSTT(J) + XSTTS(J) ) . 

This follows the same predictor-corrector sequence as the Gabutti scheme 
used in the updating routine and is thus consistent with the rest of the 
scheme. At the end of the third step, the newly corrected values of the 
shock (and connecting line) position and velocity, XS and XST, are 
substituted for the old values. 


(b) Calculation of XSN from J » JS+1 to JMAX 



(c) Blending of XSN from JWALL to JMAX 
Pig. 5.2 Construction of Line I - IS 





Chapter 6 


SCHEME 3 - FLUX VECTOR SPLITTING EULER 
6.1 Introduction 

The schemes described in Chaps. 4 and 5, when used together, will 
produce an accurate resolution of the flow field and will give the 
correct jump relations through shock waves. However, the schemes have 
several disadvantages which limit their utility. First, shock fitting 
is not a technique which is easily implemented and the need for adaptive 
gridding adds to the computer run times and complicates the coding. 
Second, since the governing equations are not solved in conservation 
form, there is no guarantee that mass or momentum is conserved through- 
out the flowfield. Third, since the equations are solved explicitly, 
there is a restriction to low maximum allowable CEL numbers and as a 
result covergence to a steady state is slow. It was this restriction 
which applies to all explicit schemes which led to the development of 
implicit schemes in the late 1970's. The CFL restriction on the Lambda 
algorithm was removed several years ago with the introduction of 
implicit Lambda schemes by Dadone and Mapolitano [43] and Dadone and 
Magi have developed a "quasi-conservative" Lambda formulation [44]. 

The scheme described in this chapter is an approximate factori- 
zation scheme which was developed by von Lavante. The scheme is 
implicit and solves the governing equations in conservation form. The 
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scheme is described in Ref. [45] as the "true Jacobian scheme" which 
refers to the implicit operator which is used. 

6.2 Algorithm Development 

The two-dimensional, unsteady, compressible Euler equations in 
conservation form are 


°t + F x + G y = 0 (6.1) 

where 

Q - Cp , pu , p v]^ 

2 T 

F = [pu, pu + p, puv] 

2 T 

6 = [p v, puv, pv + p] 

As was the case in Scheme 1, the pressure is related to the density and 
the velocities by the steady state energy equation 



a 2 - IE 

P 

These equations are nondimensionalized as follows 

x* * x/L, y' = u/L, u 1 * u/a Q , v ' 


( 6 . 2 ) 

(6.3) 





h t /a 


2 

o 
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Equations (6.2) and (6.3) then become 





(6.4) 


The algorithm to solve Eq. (6.1) first does a time discretization 
of the vector Q by a truncated Taylor series as follows 

Q n+1 • Q n + Lt + 0 (*t 2 ) (6.5) 

Let aQ = Q n+ * - Q n and express the time derivative of Q n+ * by another 
Taylor approximation and Eq. (5.5) becomes 


4Q * ^ (§^ ♦ St k (f^)) + 0 (it 2 ) 

Using Eq. (6.1), this becomes 

aQ = -At [ (Fj + Gj) + At ^ (f; + Gj)] (6.6) 

The conservation form of the Euler equations has the property that the 
flux vectors F and G are homogeneous functions of degree one of the 
vector Q. This means that 

F = AO and G - BQ (6.7) 


where A and B are the Jacobian matrices a F/a Q and aG/aQ respectively. 
These matrices are given below 


A 


0 



(y*1)u 2 ) - 


-uv 




(6.8a) 
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B 

Equation (6.6) 
form as 


1 


0 

-uv 

|(eu 2 - ( Y +i)v 2 ) 


0 

V 


1 

u 

.+1 


(6.8b) 


after some manipulation, can be expressed in operator 


i 1 + st 'It ♦ H>] ^ (H ♦ If) < 6 - 9 > 

This equation can be approximately factored as 

!■ * At ft] [I lQ . + 

= R n (6.10) 


Equation (6.10) is solved in three steps as follows. First, the right 
hand side of (6.10) is found using the flux-vector splitting as 
described in the next section. Then, aQ is found in two steps as 

follows 


r 3 A. ~ n 

[ i + it J7I 0 ■ R 

[i 4 q . n 

Finally, Q is updated. 


( 6 . 11 ) 


Q n+1 = Q n + aQ 


Generally, block tridiagonal matrix inversions must be used to solve Eq. 
(6.11). Equation (6.10) can also be transformed from cartesian to 
general coordinates 
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where 


n ♦**§ n ***3 


0 = Q/J 


■At(F + G ) 
£ n 


( 6 . 12 ) 


F = 5 F + 5 G 
x y 


G * n F + n G 
x y 


4 = x v - x y. 

0 5 n n J 5 

Since the Euler equations are hyperbolic, the Jacobian matrices 1 
and B can be diagonalized by the following similarity transformations 

.-1 

(6.13) 


A = M A. M 
A 


B = N A d N 

b 


-1 


where M and and N are the right eigenvector matrices and a a and a b are 
the diagonal matrices of the eigenvalues. It can be shown by examining 
the equations in primitive variable form that the eigenvalues of A are 


X2 = u(y+1)/2y + SQ (6.14) 

X 3 = u(y+1)/2y - SQ 

where 

SQ = [(u(y-1)/2y) 2 + a 2 / Y ] 1/2 


In Refs. 46 and 47, von Lavante describes how the Jacobian matrices 
can be diagonalized to permit even greater computational simplification. 
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6.3 Flux^Vector Splitting 

The flux vectors F and G are each split into forward ( + ) and 
backward (-) fluxes by a method introduced by van Leer [481. To use the 
flux vector F as an example (the requirements for G are the same), the 
flux F is split as 

F (0) - F + (Q) + F“(0) (6.15) 

The second requirement is that 

A + = 3 F + /3 Q have all eigenvalues > 0 

A" = 3 F”/ 3 Q have all eigenvalues < 0 

The flux-split components F + and F - must be continuous and satisfy 

F + = F for Mach No. fl > 1 

F“ = F for M < -1. 

The components are further required to correctly model the symmetry of F 
with respect to the Mach No. M such that 

F + (M) * t F'(-M) if F(M) « ± F(-M) 

The split Jacobian matrices A + and A” must be continuous at sonic and 
stagnation points. This requirement is important as other types of 

splittings do not accomplish this and these splittings produce 
oscillations when the eigenvaluve change signs. Next, it is required 
that for subsonic flow, A + and A - must each have one eigenvalue 
vanish. Thjs requirement makes it possible to capture shocks with no 
more than two interior cells. Van Leer satisfies these requirements by 
appropriate choices of polynomials to represent F(M) and G(M). 
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6.4 Spatial Discretization 

The right hand side of Eq. (6.12) can be integrated over a finite 
volume which consists of the interior of the quadrilateral cell with the 
grid points at (i, j), (i-1, j), (i, j-1), and (i-1, j-1). Thus, using 
Green's Theorem, we have 

aQ = -At / (F + G ) dA 
e A 5 n 

= -At J F dn - % dc (6.16) 


The line integral above can be approximated as 



-At 


‘V 


1/2 " F i-l,j-l/2 


where, for example, F. . , 

i , j-1/2 

with midpoint (i, j-1/2). Since the flux F is split into F + and F“, 
each component must be evaluated at the appropriate cell face. The 
method used to evaluate the fluxes is the MUSCL type differencing of van 
Leer [49]. In this approach, the fluxes are extrapolated to the cell 
faces according to the signs of the eigenvalues. Thus 

n - < 3 F !- f m > /2 


Vl/2,j " Vl/2,J-l) (6-17) 


is the flux of F through the cell face 


F i+l/2 = F i+1 ' F i+2 )/2 


(6.18) 
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The values of 0, F and G in cell — ( i , j ) are considered as representative 
cell averages. 

In the implicit operator (i.e the left hand side) of Eq. (6.12), 
the Jacobian matrices 1 and B are also split such that 

~ + 

A = a F /a Q 
A" = 3 F""/3 Q 


Equation (6.12) thus becomes 

[ I + At(3 + 3 V)] [I + At(3 'B + + 3 TT)1 aQ = RHS (6.19) 
5 £ n n 

The derivatives in Eq. (6.19) are taken as one-sided forward or backward 
differences depending upon the sign. Thus for example 

[I + At (a B + + 3 B")] aQ * aQ. + At(A + aQ. - A. + aO. .) 
n n 1 ii i-l i-l 

+ At(A. +1 a0. +1 - A. aQ.) (6.20) 

The form of Eq. (6.20) results in block tridiagonal matrix inversions to 
solve. As the steady state is approached, both sides approach zero, the 
first order accuracy on the left hand side of (6.18) does not affect the 
accuracy of the right hand side which is second order. 
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Chapter 7 


SCHEME 4 - UPWIND NAVIER-STOKES 
7.1 Introduction 

The full compressible Navi er- Stokes equations were solved using an 
upwind approximate factorization scheme for test problems 1 and 2. The 
code developed by Rumsey of NASA Langley, is fully vectorized to run on 
the CDC VPS 32 supercomputer, and is accurate for unsteady flows [50]. 
Both the laminar and Reynolds averaged turbulent Navier-Stokes equations 
were solved. A Baldwin-Lomax turbulence model was used and the Reynolds 
number was varied from 10,000 to 100,000. 

7.2 Governing Equations 

The full set of conservation equations in two-dinensions includes 
conservation of mass, conservation of momentum and conservation of 
energy. The momentum equations are the compressible Navier-Stokes 
equations as derived In Chap. 3. In vector form, the equations are 

VV® = ir<V~ s > o.n 

t £ n R g C n 
where Q, F, G, R and ? are 

Q = Q/J 

F = (e x F + G)/J 
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G = (n + n G)/J 
x y 

~R = (c x R + C y S)/J 

3 * (n R + n S)/J 
x y 


and 0, F and G are given by 

Q = [p , pu, P v, e] T 


F = 


pu 

2 

pu + p 
p uv 

{e + p) u 


- kT. 


G = 


pv 

puv 

2 

p/ + p 

(e + p) v - kT 


R and S are given below as 

R 


T xx’ T xy* ^4^ 


s = CO, T xy , T yy , S 4 ] 


where the components of R and S are 


xx 


(X + 2b) (5„ u ? ♦ u n ) + x( 5y + n y »„) 


yy 


(X * 2„) (5 y * { + * y »„> *X(5 X u { + „ x x n > 


= n (e u+n u+e v+n v ) 
xy wv *y e y n *x 5 n x n 


R 4 ' * u T xx * VT xy + “ (? x \ ^t/CPrly-l)] 


u T xy + VT yy * " ( 'y T 5 + "y T n )/[Pr(l, ' 1)3 



86 


2 

The Stokes hypothesis that — x * - -j u is assumed and the Prandtl 
numbers is taken as a constant. 


7.3 Numerical Algorithm 


Equation (7.1) is solved using upwind approximate factorization as 
described in references 29 and 50. The algorithm is given as 


[I + At(3~ + 3 + ^") - R _1 At 3 (M/J)l 

55 e 5 

x [ I + A t (3 ” 1* + 3 + B‘) - R -1 At(3 (N/J) 1 A^ n 
n n e n 

= - At R n 

n • 'v'f ^ *** • ■ + #v- 

R = 3 F+3 F + 3 u + 3 G 

5 5 n n 

- R _1 (3 *R + 3 'S) 
e 5 n 


(7.2) 


(7.3) 


and where as before 


_ aF 1 

‘ 3^ 

. 3G* 

‘ 5TT 


A and P are the Jacobian matrices given as Eqs. ( 6 . 8 a) and ( 6 . 8 b) and 


'h =5 A + 5 B 
x y 

B - q A + n B 
x y 

The (+) and (-) superscripts on the F and 6 terms indicates the flux 
splitting which is done according to van Leer [48]. The (+) and (-) 
supescripts on the partial derivatives terms denotes the direction of 
differentiation. Thus, for example, d~ is a backward difference. 
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while 3* is a forward difference. All viscous terns are centrally 
differenced. The M and N matrices are linearized viscous terms and 
are given by Steger in reference [51]. The scheme is first order 
accurate in tine and second order accurate in space. 

7.4 Turbulence Model 

The turbulence model that is used in this scheme is a two layer 
algebraic eddy viscosity model developed by Baldwin and Lomax [52]. The 
model computes an eddy viscosity y t which is then added to the 
molecular viscosity y to get the total viscosity. The model follows 
from a previous model developed by Cebeci [53] but avoids the necessity 
for finding the edge of the boundary layer. It has been shown to give 
good results in separated flows and in wake regions. 

In the inner layer y^ is computed as 



The length i is obtained using the van Driest formulation as follows 

l = ic y [1 - exp(- y + /A + )] (7.5) 

< is the von Karman constant equal to .4, y + is the non-dimensional 
wall unit 

y + = y(T y ^)* 5 /y (7.6) 

and A + is equal to 26. 

In the outer region, the turbulent viscosity is given by 
w t ” * ^cp * ^wake * ^Kleb ^ 


(7.7) 
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where 


and y nax 


and 


K * Clauser constant = 0.0168 


C cp ■ 1.6 


^wake ” f1l ' n ^ Y nax p max* C wk Y max ^max^max-^ 


r K]eb (y) ■ Ci + s.b (c„ loh /Y mav n 


6-.-1 


Kleb' max J 


and F max are the location and maximum functional value given by 


F(y) = y (u) [1 - exp(-y + /A + )] 


C wk - °‘ 25 
C Kleb = 0,3 * 



Chapter 8 


INITIAL AND BOUNDARY CONDITIONS 
8.1 Introduction 

As stated earlier, the equations of fluid dynamics pose an initial - 
boundary value problem which are then solved numerically. Initial 
conditions (i.e. the initial state of the fluid) must be supplied for 
the solution to proceed. It is usually assumed that the final steady 
state solution is independent of the initial conditions, although this 
is not necessarily the case and the author is unaware of any 
mathematical theorems that state this is true. 

On the other hand, the boundary conditions are crucial to achieving 
a correct solution of the problem and must be consistent with the 
physics of the problem. The boundaries of the computational domain 
include both physical and artificial boundaries. Physical boundaries 
are normally walls and it is generally, although not always, possible to 
specify conditions at these boundaries in a straightforward manner. 
Artificial boundaries exist due to the necessity of having a finite 
computational domain and the specification of conditions at these 
boundaries is more difficult and open to question. 

The boundary conditions are also important insofar as the stability 
and convergence, properties of the numerical scheme. Improper specifica- 
tion of the boundary conditions can lead to instability or slow conver- 
gence of the computations. 
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8.2 Initial Conditions 

The Initial conditions used in each of the four problems varied 
depending on the scheme used and are described in the following 
sections. 

8.2.1 Problem One 

The initial condition used in scheme one was obtained from the 
complex velocity potential. It is shown in texts of ideal fluid flow 
that the potential and stream function together form an analytic 
function which is called the complex potential as follows 

F(c) - 0 + i * 

- A c (z) (8.1) 

where A is a constant. The derivative of the complex potential is then 

F'(c) = dF/dz 

= u - iv (8.2) 

The right-hand side of Eq. (8.2) is called the complex velocity and is 
the complex conjugate of the velocity vector u + iv. For problem one, 
the complex velocity is given by 

F'(c) - Aw C( c - 1)/U + 1)] ’ 5 

= ^ = u + iv (8.3) 

In order to have u go to u^ as c goes to infinity, it is necessary to 


have 
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A~= u /it 
00 

Equation (8.3) then becomes 

u + iv « u [(c - l)/(c + 1)] * 5 
00 


(8.4) 


The value of u^ in this case is specified according to the pressure at 

the exit of the channel since there is no "infinity." Isentropic flow 

is assumed and the speed of sound and the Mach No. at the exit are 

obtained from the exit pressure. The velocity at the exit is then taken 

to be u . The pressure is obtained by assuming the flow to be 
00 

isenthalpic and using the initial u and v components in Eq. (4.3) to get 
the speed of sound which is then related to the pressure by the 
isentropic relation 


p = Ca 2 /y] 


1/(1-y) 


The initial condition used in scheme two is a converged solution 
from a previous run of scheme one. First, the solution is scanned along 
lines of constant n to determine the shock location. The shock is 
initially placed at the midpoint between the upstream supersonic and 
downstream subsonic points. In order to avoid kinks in the initial 
shape of the shock, which result when the e coordinate value of the 
midpoint shifts, the shock shape is then smoothed by a simple iterative 
routine which keeps the shock between the supersonic and subsonic 
points. 

Once the initial shock shape has been determined, a new shock 
fitted grid is introduced and the old solution is interpolated onto the 
new grid. The distribution of points along the lines of constant n is 
done by the following second degree polynomial: 
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e = e o + ( £ 1 * E o^ ^ a l * + 
where aj and a 2 are calculated to satisfy 


X 

II 

o 

at 

e 

= e o 

X = IM-1 

at 

e 

= E 1 

X = IS 

at 

e 

= 


IS is the integer value of the shock point and e $ is the e location 
of the shock. 

The interpolation from the old grid onto the new one is done by a 
first order linear interpolation. The values of the flowfield u, v, and 
P at the upstream shock point are determined by a linear extrapolation 
from the two previous upstream points on the old grid. The Rankine- 
Hugoniot jump relations are used to obtain u, v, P and the entropy S on 
the downstream shock points. Since scheme two uses the entropy as the 
energy variable and since this is not one of the variables used in 
scheme one, the entropy at each point is obtained from the pressure and 
the speed of sound using Eqs. (4.4) and (5.3). The normalization for 
the entropy sets the entropy of the incoming flow equal to zero and so a 
calculated value of the entropy greater than zero at a point indicates 
that the fluid has passed through a shock or some other dissipative 
mechanism. A calculated entropy of less than zero is erroneous and is 
instead set to zero. 

The initial conditions used in scheme three are obtained from the 
specified Mach number at the entrance ancf exit of the channel and the 
assumption that the v component of velocity is zero at these locations. 
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Then using the assumption that-the flow Is isenthalpic the speed of 
sound, the pressure, density and the u component of velocity is obtained 
at the entrance and exit. The points between the entrance (1=1) and 
the exit (I = IM) are then determined by a linear interpolation based on 
the I coordinate. 

The Navier-Stokes code (scheme four) was used in a "flate plate" 
option which assumed the upper computational boundary is a free stream 
boundary instead of as a wall as the first three codes assumed. The 
code, therefore, initialized all flow quantities to the specified free 
stream values. 

8.2.2 Problem Two 

In scheme one, the complex velocity is used as in Sec. 8.2.1 and it 
is found that the velocity is given by [29] 

u + iv = u (1 - 1/z 2 ) (8.6) 

00 

and the pressure is given by Eq. (8.5) 

The initial conditions in scheme two were obtained using the same 
procedure as in problem one. 

The initial conditions in scheme three were obtained using the same 
procedure as in problem one. However, since the inflow and outflow Mach 
numbers are the same, the initial condition for this problem is one of 
uniform flow. 
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8.2.3 Problem Three — 

For scheme one, it can be shown that the velocity is given by 

u + iv = u /[n(dz/d;)] (8.7) 

00 

where dz/dc is given by 

dz/dc = n(z 2 - l)/( c 2 - 4) 

and n is given by Eq. (2.8). Again, the pressure is found from the 
assumption that the flow is isenthalpic and isentropic as in problems 
one and two. 

The initial conditions used in schemes two and three are the same 
as those used in problems one and two. 

8.2.4 Problem Four 

The initial condition used for the airfoil problem was developed by 

(1) first setting surface where the velocity was set equal to the 

component of the uniform free stream u tangent to the wall and then, 

00 

(2) using a simple laplacian operator to smooth the flowfield prior to 
beginning with the actual flow solver. The reason for using this 
procedure was to minimize transients which could cause instability in 
the initial solution and which would be due to initial sharp gradients 
that would exist such as at the leading edge where the velocity would be 
zero at the wall and free stream one point off the wall. 
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8.3 Wall Boundary Conditions 

The condition of no flow through the wall was applied at all wall 
boundaries. For the inviscid calculations, this is the basic condition 
whereas for the viscous calculations, the additional condition of zero 
velocity tangent to the wall or the so-called "no slip" boundary 
condition must also be imposed. The following sections described in 
detail how these conditions were implemented in each scheme. 

8.3.1 Scheme One 

Referring to Eqs. (4.21), the no penetration boundary condition 
implies that V = = V t = 0. The two governing equations reduce to 

P + + .5 [x + P* + a' p; + fi + P + + ft" P* 

t 5 £ n n 


— (A + U + 
iV J 5 

- a" 

u; + n + V + - n" V)] = P c 
C n n s 

(8.8a) 

.5 [\ + U + + A" 

If + 

(A + P + - A - P")] = U e 

(8.8b) 

5 

5 

Y n n s 

.5 [n + v + + n“ 
n 

V" + 
n 

— (n + p + - ft" p‘)] = v c 

Y n n s 

(8.8c) 


Equation (8.8b) can be used as it is, but before Eq. (8.8a) can be used, 

it must be modified to eliminate the terms ft + P + and ft + V + at the 

n n 

lower wall and the terms n' P" and ft" V" at the upper wall since 

n n 

they require information from outside the computational domain. This is 
done by combining Eqs. (8.8a) and (8.8c) to get 



r 
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p f + .5 [x + P* + x‘ p; + 2 {r p + 

t 5 5 n 


+ — (x + u + - x" u‘ + 2 n + v + )] 

a/T* 5 5 n 


+ — 1 — 
a/ J 


V 


s 


which applies at the upper wall and 

P t + .5 [x + P* + x“ P" + 2 n‘ P" 
1 5 5 n 


(x + u + - x" U + - X* u: - aT V’)] 
5 5 5 n 


p s - 


a/ J 4 4 4 " a/ J 

which applies at the lower wall. The governing equations can therefore 
be used at the walls and the three explicit scheme used as in the 
interior points. 


8.3.2 Scheme Two 

The same concept as just described for scheme one is also applied 
in this case. The resulting equations are slightly different and will 
not be given. 

8.3.3 Scheme Three 

Since scheme three solves the governing equations implicitly, the 
implementation of the scheme at the boundaries poses difficulties and 
explicit boundary conditions are used instead. The velocity V and 
density at the wall are first found by extrapolation from the inside. 
The pressure at the wall is next found by extrapolation from the 
interior using the momentum equation in the n direction 

§£ ■ p (u Z ♦ v Z )/[R(n^ 4 ^)- 6 ] 


(8.9) 
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where R is the radius of curvature and is given by 

R = (1 + {f*) 2 ) 1,5 /f" 

f' = dy/dx (8.10) 

The density at the wall is now recalculated by using the assumption that 

p • y P Y (8.11) 

which is not valid if there are entropy gradients normal to the wall but 

which nay be regarded as approximately true. Using this assumption, the 
density gradient along n Q is given by 

f^<l/. 2 )f£ (8.12) 


8.3.4 Scheme Four 

As stated earlier, for the viscous calculations, both the u and v 
components of velocity are taken to be zero to satisfy the no-slip 
boundary condition. The density at the wall is next obtained by a zero 
order extrapolation from the point above it. Next, the speed of sound 
one point off the wall is found by using the energy equation as follows. 
The density at the same location is then used in conjunction with the 
speed of sound to get the pressure at the wall. The pressure at the 
wall divided by (y - 1) is then the total energy which completes the 
procedure. 
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8.4 Subson+c Inflow Boundary 

This section and the following two involve "artificial" boundaries 
in the sense that they exist in the computational domain but not in the 
physical domain. Since boundary conditions must be imposed to make the 
problem solvable and since these conditions must reflect the actual 
physical situation, the correct specification of these conditions is 
fairly important. Much of what has been done up until the present time 
involves making some simplifying assumptions and frequently 
characteristic theory is used. At the upstream boundary, the flow is 
subsonic and entering the computational domain. If the flow is assumed 
to be inviscid so that the Fuler equations apply, the governing 
equations in primitive variable form are 


Q t + A Q x + B Q y = ° 


where 



(8.13) 


(8.14a) 


A = 


u 


0 

0 

0 


P 

u 


0 

YP 


0 0 * 

0 1/p 

u 0 

0 u 


(8.14b) 


B = 


v 

0 

0 

0 


0 p 

v 0 
0 v 
0 Y P 


0 

0 


1/p 


v 


(8.14c) 
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The eigenvalues of A are u, u, a- + a and u - a. The eigenvalues of B 
are v, v, v + a and v - a. The directions u + a and v + a are referred 
to as bicharacterisitics and are the intersections of the characteristic 
monge conge with the planes x = constant and y = constant. 

It is well known that for purely one-dimensional unsteady isen- 
tropic flow, the equations of gas dynamics possess Riemann invariants 
which are constant along characteristics. These functions are 

R + = u + 2 a/(y - 1) along dx/dt = u + a (8.15a) 

R" = u - 2 a / (y - 1) along dx/dt = u - a (8.15b) 

At a subsonic inflow boundary, the R + value is constant along the 
characteristic entering from outside of the computational domain and can 
therefore be specified and used in conjunction with R" from inside the 
domain to get new values of u and a. Equation (4.13) does not possess 
Riemann invariants but according to the theory of Kreiss, the number of 
conditions which must be specified at an inflow (or an outflow) boundary 
in order to have a numerically stable boundary condition must agree with 
the number of characteristic lines that approach the boundary from the 
outside [52]. If the flow is assumed to be at least locally one- 
dimensional, this would mean that there must be three quantities 
specified from the outside and one extrapolated from the inside at the 
inflow boundary and the reverse at the outflow boundary. 

8.4.1 Schemes One and Two 

The conditions which were assumed for problems one through three 
were that the contravarient component of velocity V is zero and that the 
flow is isentropic and isenthalpic at the inflow boundary. The first 
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assumption is an approximation to- the actual physical situation and is 
not exactly correct while the thermodynamic assumptions are correct in 
the steady state. Since V and V t are assumed zero, the governing 
equations are 

p t * .5 [x + p; ♦ X- P- ♦ n* p; ♦ n- p; ♦ -X- <x + U* - X- up] • P s 

dy u 

U t * 5 C>+ U x n ’“i* (x+ P X • x' Pp] =U S 

“V t 4i T ( “‘ P y ‘ n ' P y'> ■ v s 

Since the flow is subsonic, it is not possible to use these equations 
since P and U require information from outside the domain. But the two 
equations for P t and U t can be combined to eliminate these terms to get 

P t - — u t + x" p' + .5 (n + P* + o' pp + —I— (u . X- up = P c 

‘ ‘ x (8.1*) y •rr s 

It is possible to relate P t to U t by the 1-D energy equation 

a 2 /(y - 1) + .5 (u 2 + v 2 ) = h Q = const (8.17) 

It is not necessary but if it is also assumed that v = 0, then 

a ■ U/J* 5 (8.18) 

and since the flow is isentropic 


a 2 = y P^- 1)/T 



101 


which leads to the equation — 

P t - - y U U t /(a 2 J) (8.19) 

which can then be substituted back into Eq. (8.18) to get a single 
equation for U t as follows 

„ u s ->~ u;-a J- 5 [P s -x- p;-.5te* p;*a-p y )] 

t [1 + u/a J ,S ] 


P t can then be obtained using Eq. (8.19). 

8.4.2 Scheme Three 

The inflow boundary condition used in this scheme is described by 
von Lavante in [47] and will only be briefly discussed here. The 
assumption is made that the flow is locally one-dimensional so that the 
governing equations can be written in the characteristic form 

Hr + A f£ * 0 < 8 - 20 > 

« 

where W is the characteristic variable given as 



V 


Wj ■ 1 + u(y - l)/2 so 1D 

W 2 = 1 - u(y - l)/2 SQ 1D 


A = 



0 

x 2 


( 8 . 21 ) 


„ SQ X - C(u( Y - 1 ) / 2 Y ) 2 + a 2 /Y]* 5 
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The characteristic variable propagates along the characteristic dx/dt 
= u + a and is therefore specified whereas the variable W 2 propagates 
along the backward characteristic dx/dt = u - a and is taken from the 
inside of the computational domain. Equations (8.21) are then solved 
simultaneously on the boundary to get the Mach number at the inflow 
which is 

M 2 = l/{ (y - l) 2 [1/(W 1 - W 2 ) 2 - 1/4]) (8.22) 

The Mach number is then used to calculate p, u, v and p from the 
isentropic relationships and the energy Eq. (8.17). The procedure is 
very similar, at least in principle, to a "time split" inflow boundary 
condition Introduced by Tong [55] and avoids the necessity for assuming 
that v = 0 at the inflow boundary as was the case in schemes one and 
two. Furthermore, it has been found that the procedure gives good 
results even when the upstream boundary is close to a leading edge. 

8.4.3 Scheme Four 

The condition used in this situation first assumes the v component 
is zero so the flow is taken to be locally one-dimensional. The 
pressure at the inflow boundary is then found from the energy equation 
using the known values of the total energy and the u and v components of 
velocity at the next point downstream. The pressure is then used to get 
the density assuming the flow is also isentropic. Finally, the u 
component of velocity is found using the isenthalpic energy Eq. (8.17). 
This procedure should thus allow pressure waves which are traveling 
upstream,to escape without reflection. However, the assumption that v = 
0 at the boundary is not strictly correct and may introduce errors. 
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8.5 Outflow Boundary Condition 

At the outflow boundary, there are three characteristics from 
inside the computational domain and one from outside and hence according 
to the theory of Kreiss, one condition may be specified at the outflow 
boundary. This condition is usually taken to be the pressure since this 
determine the flow through the channel at the steady state. For 
example, in a De Laval nozzel, the backpressure completely determines 
the flow inside the nozzel and the conditions at the throat. However, 
if the pressure is specified as constant at the outflow boundary, then 
during the convergence from the initial condition to the final steady 
state solution, pressure waves which may develop inside the 
computational domain are reflected at the outflow boundary instead of 
passing through and this will slow convergence. 

The approach used in the first three schemes was introduced by Rudy 
and Strikwerda [54] and was based on earlier work by Enquist and Majda 
[55] and Hedstrom on nonreflective boundary conditions. The idea is to 
apply the following equation at the outflow boundary 

ft “ p c fr + a {p " p e ) (8 ' 23) 

where a is some constant and p e is the specified pressure at the 
boundary. T « the steady state, the pressure p should be equal to p e . 
The numerical approximation to Eq. (8.23) is 

p n+1 * [p n + a At P e + p n a 11 (u n+1 - u n )]/(l + a At) 

where u n+ * is extrapolated from the inside to the boundary. This 
procedure has been found to work well but has the disadvantage that the 
parameter a must be chosen. As a approaches zero, the influence of the 
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specified back pressure is also relaxed and approaches zero. As a is 
increased, the outflow boundary becomes increasingly reflective and p is 
maintained close to the specified p e . Rudy and Strikwerda recommend 
that a be chosen to optimize the convergence to the steady and show 
that a has a significant effect on the number of iterations required 
for a ty r cal test case to converge. Figure 8.1 shows the results of 
applying this boundary condition to problem one. 

In scheme four, the pressure was simply held constant at the 
outflow boundary while the density and the u and v components of 
velocity were set equal to the next upstream point. This is an entirely 
reflective boundary condition but was used in the flat plate option 
because the flow assumed to be mainly all boundary layer and hence 
parabolic. This assumption was not true in the cases modeled but is the 
physical condition imposed is correct in the steady state. 

8.6 Far Field Boundary Condition 

The far field refers to the flow field at a distance away from the 
body which is significantly greater than the length scale of the body. 
For example, it is common practice to put the outer boundary for two- 
dimensional transonic calculations about airfoils at least ten chord 
lengths away from the airfoil. The reason for this is that one would 
like to specify conditions in the free stream where the flow is 
unaffected by the body but this only occurs at an infinite distance away 
from the body. Since it is impossible to construct grids which go out 
to infinity, it is necessary to specify boundary conditions in the far 
field where there are still perturbations about the free stream 


conditions. 
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In problems one and two, the upper boundary was treated as a solid 
wall in all of the schemes except scheme four. Since this scheme solves 
the viscous equations, it was deemed desirable to avoid treating the 
upper boundary as a wall since doing so would necessitate grid 
clustering at the top and lower boundaries due to the existence of 
boundary layer growth on both surfaces. Therefore, this boundary was 
treated in exactly the same manner as the outflow boundary was treated. 
The pressure was held constant at the free stream pressure and the 
remaining variables were extrapolated from inside the computational 
domain. This approach may not produce the best convergence, but the 
results appear reasonable at the top boundary. Attempts were made to 
develop characteristic based nonreflective boundary conditions for the 
top boundary for schemes one, two and three which would permit this 
boundary to be relatively close to the lower boundary and allow for a 
better comparison between these schemes and scheme four but these 
efforts were not successful. Instead, the upper boundary was 
progressively moved away from the lower wall until the perturbations at 
the upper boundary were reduced to an acceptable level. 

For problem four, the outer boundary condition which was used was 
developed by Thomas and Salas [56]. The approach assumes that wavelike 
transients in the physical domain arrive at the outer boundary in a 
direction mainly normal to the boundary as shown in Fig. 8.2. The 
flowfield at the outer boundary can be thought to consist of a uniform 
flow plus small perturbations due to these transient waves which in time 
should decay to the final deviations from the uniform flow. Therefore, 
the gradients along the outer boundary are assumed small in comparison 
to the gradients normal to the boundary the following characteristic 
relations are presumed valid 





Vlavelike transient 
approaching boundary 


8.2 Far Field Boundary Condition. 
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dS/dt = 0 

■along dn/dt = U 

(8.21a) 

dV/dt = 0 

along dn/dt = U 

(8.21b) 

dU 1 d£ _ n 
3T + • 0 

along dn/dt = U + a J*^ 

(8.21c) 


where S is the entropy and U and V are the contravarient velocities in 
the 5 and n directions respectively. Since the flowfield at the 

outer boundary Is isentropic for inviscid flow and nearly so even for 
viscous flow, Eq. (8.24) implies the existence of Rienann invariants at 
the outer boundary and the equivalent of Eq. (8.15) in a direction 
normal to the outer boundary is 

R" = • n + 2 a/( Y - 1) 

along dn/dt * ^ • n + a (8.25) 

where 

1/ . n = V/tJ (x* ♦ y*f 5 ] 

The invariant R", which reaches the boundary from outside the computa- 
tional domain, Is specified from the free stream values of u and a 
applied at the outer boundary and the values of R + is taken from inside 
the computational domain. The equations for R + and R" are then solved 
simultaneously to get the new values of V and a on the boundary. 
Depending on the sign of V, the values of the remaining variables U, p, 
and p are then either taken from inside or outside the computational 
domain. If V is negative, the flow at the, boundary is in and the free 
stream values of u, p, and p are used. If V is positive, the flow at 
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the boundary is out and U, p, “and p are extrapolated from the inside 
of the domain. This procedure has been found to give good results in 
the computations performed as a part of this investigation. However, 
Roe has questioned the practice of assuming locally one-dimensional flow 
at remote boundaries and claims that it does not result in the monotonic 
decay of the pressure and the radial component to the expected free 
stream values [57]. More work will undoubtedly be done on the treatment 
of far field boundaries in the future in view of the importance of the 
topic. 



Chapter 9 


RESULTS AND DISCUSSION 
9.1 Problem One 


9.1.1 Inviscid Results 

The three inviscid codes were run using various values of n for the 
lower and upper boundaries, referred to as n 0 and respectively. The 
lower boundary n Q was varied to study the effect of changing the 
geometry from a channel with a gradual one-sided expansion to a channel 
with a one-sided sudden expansion (rearward facing step). It was found 
that a value n Q = 2 produced a shock at the expansion corner near 
5 = -1 and the flow was attached. As n Q was reduced, the flow would 
eventually separate (usually before n Q = 1) and as n 0 was further 
reduced the recirculation zone would grow in size. The effect of moving 
the upper boundary was also studied; it was found early in the 
investigation that the position of n^ strongly affected the flow in the 
entrance region of the channel ahead of the expansion corner. This is 
because the ratio of the entrance area to the exit area is increased and 
approaches 1.0 as nj is increased and the result can be seen in Fig. 9.1 
which shows the Mach number distribution across the inflow boundary 
for nj = 20, 40, and 100 using scheme one. As is increased, the 
inflow Mach numbers are reduced as would be expected. 
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The three inviscid codes “were run on 31x41 grids of which an 
example (with n Q s 2 and = 20) is shown in Fig. 9.2. In addition, 
the shock fitting was also repeated using 161x81 grids to determine the 
effect of refining the grid by a factor of two. In general, scheme one 
was run for 3000 iterations in each case, which usually resulted in a 
residual drop, measured by the nom 3 t p » between two and three 
orders of magnitude. The slow rate of convergence is due to the CFL < 1 
limitation and the fact that no convergence acceleration methods other 
than local time stepping were used. 

The shock fitting runs were done in continuous sequences with n 0 
initially equal to 2 and then gradually reducing n Q in a step-like 
manner. Thus, the lower wall was first brought to n 0 of 1.8 and this 
position held constant until the shock movement was reduced to an 
acceptable level. This procedure was repeated by reducing n Q in steps 
of 0.2 and usually at least 2000 to 3000 iterations were required on the 
coarse grid to steady the shock, while on the fine grid 4000 to 5000 
iterations were required. 

The implicit flux-vector splitting scheme was run for 2000 
iterations in most cases and the residual drop, measured by the norm of 
a t p, was usually at least three orders of magnitude. In some cases, a 
constant global time step was used, while in other cases a constant 
local maximum CFL generally between three and six was used. It was 
found that the later approach gave the best convergence results and was 
therefore used for most of the later runs. 

The Mach number on the upstream side of the shock foot is shown in 
Figs. 9.3, 9.4 and 9.5 for positions of 20, 40 and 100 respectively. 
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In Fig. 9.3, it is apparent that scheme one produces significantly 
higher Mach numbers as n Q is reduced than the other two schemes. This 
is believed due to the nonconservative formulation of scheme one. There 
is fairly good agreement between the coarse and the fine grid solutions 
obtained using scheme two. The Mach numbers obtained using the flux 
vector splitting scheme were consistently lower than those obtained 
using schemes one and two, perhaps due to greater dissipation at the 
shock in scheme three. The lowest n 0 for which solutions were obtained 
was 0.4 for schemes one and two and 0.1 for scheme three. It was only 
possible to get a solution at n 0 = 0.1 for scheme three by changing the 
grid in steps to gradually reduce n 0 » while at the same time increasing 
cluster points at the expansion corner. It can be seen that the Mach 
number at the shock foot reaches a maximum for schemes two and three at 
n Q = 1.2, and as n 0 is further reduced, the Mach number gradually 
becomes lower and the results using scheme three suggest it would go to 
a 1 imit of 1. at n Q = 0. 

The results at = 40 are shown in Fig. 9.4. The results using 
scheme one, again, show that as n fl is reduced, the Mach number at the 
shock foot increases. The results from the shock fitting scheme two 
show two different tendencies on the coarse and fine grids. On the 
coarse grid, the trend is similar to scheme one, whereas on the fine 
grid the results are similar to those of scheme three. The anomalous 
results using scheme one were checked by doing this series of runs over 
again with 4000 iterations being taken when n Q is reduced each 0.2. The 
results were well converged and do not appear to be in error. The 
increase in the Mach number is partially explained by the shock position 
as will be discussed later. Again, the results using scheme three 
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show the lowest Mach numbers “and are similar to those obtained at 
Oj = 20. The lowest n 0 for which a solution was obtained was 0.1 using 
scheme three and this was again achieved by changing the grid to 

progressively cluster the grid about the expansion corner. The trends 
show the Mach number ahead of the shock decreasing to a limit of one as 
n 0 approaches the limit zero which suggests that the shock would 
eventually vanish. The results obtained using scheme three for 
= 40 agree closely with the results obtained with the same scheme for 

nj = 20. 

From Fig. 9.5, it is seen that the Mach number at the shock is 

always less with n 1 = 100 than each corresponding case at rij = 40. 

This is due to the lower velocities in the flow field ahead of the 

shock. The fine grid shock fitting results show the Mach number to be 

gradually increasing; whereas, for the flux vector splitting, the 

maximum value of M occurs at 1.2 and then decreases. No reliable 

results were obtained for n values lower than 0.4 because the resolu- 

o 

tion of the grids at the corner at {■ = -1 was inadequate. 

Figures 9.6, 9.7 and 9.8 show the shock location $ s in terms of 
the 5 coordinate along the lower wall. A value of 5 equal to -1 
corresponds to the location of the upper corner as shown in Fig. 2.1. 
The results obtained using scheme 1 are identical at * 20 and 

n j = 40 and are not appreciably different at n 1 = 100. The shock 

fitting results and the flux vector splitting results at = 20 show 
general agreement with 5^ = 40 on the coarse grid show the shock moving 
progressively towards % = -.5 as n Q is reduced. This partially explains 
the Mach numbers at the shock the position 5 = -.5 is considerably 
further around the expansion corner than the fine arid results which 
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means that the flow has accelerated more and is moving at higher 
velocities. The flux vector splitting results at all three values of 
nj agree fairly closely. In particular, it is noteworthy that the shock 
position at n Q = *1 is approaching g * -1. The shock fitting scheme 
two for nj = 100 also shows g going to a limit of -1 as n Q is reduced. 
It thus appears that this is the proper limiting position of the shock, 
but this is not achieved in all of the case computed. 

The maximum entropy in the flow on the downstream side of the shock 
is shown in Fig. 9.9 and the vorticity in the flow downstream of the 
shock is shown in Fig. 9.10 both for cases when n 1 = 20. It can be 
seen that there is a fairly consistent trend for the entropy to increase 
and the vorticity to decrease toward larger negative values as n 0 is 
reduced. This result must be interpreted with caution. As stated 
earlier, it may be the vorticity in the initial unsteady transient flow 
which leads to separation. Theoretically, none of the entropy and 
vorticity produce in the flow in the steady state actually enters the 
separation zones, when they occur. Nevertheless, it is indicative of 
the conditions which may have led to separation. 

In Fig. 9.11, the minimum value of the Mach number multiplied by 
the direction of the flow along the lower wall is plotted against for 
rij = 20. This indicates the approximate value of n 0 at which the flow 
separates for each scheme. It can be seen that scheme one produces 
separation at a higher value of (n Q = 1.68) than scheme two 
(n Q = 1.52) and that scheme three requires n 0 equal to approximately 
1.12 before separation occurs. Thus scheme three appears to be the 
least susceptible to separation as compared with schemes one and two. 
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The reattachement 5 locations- for the same set of cases is shown in 
Fig. 9.12. There appears to be general agreement in that as n Q is 
reduced, the separation bubbles grow in length and reattach further 
downstream. Similar results for the case when = 100 are shown in 
Figs. 9.13 and 9.14. It can be seen that, in general, having the upper 
boundary further away reduces the n 0 at which the lower wall must be 
before separation occurs and the 5 at which reattachment occurs is less, 
which means that the separation bubbles are not as long. 

A comparison of the result using schemes one, two and three with 
nj * 20 and n Q = 2.0 is shown in Figs. 9.15 and 9.16. Figures 9.15a- 
9.15c show Mach contours and it can be seen that the results are all 
similar. The Mach numbers multiplied by the sign of U along the lower 
wall are shown in Fig. 9.16 and again the close agreement is evident. 
The shock location is nearly identical in all three cases while the 
maximum computed Mach number reached at the upstream side of the shock 
is higher using shemes one and two than three. A similar comparison for 
the case with nj = 20 and n Q * 1.2 is shown in Figs. 9.17 and 9.18. 
Again general agreement is evident although the Mach contours of the 
flux vector splitting results appear smoother. A comparison of the Mach 
number multiplied by the sign of U is shown in Fig. 9.18. The flux 
vector splitting results, while close to separation, have not separated 
whereas the results of both scheme one and two have. It is also 
apparent that the jump in velocity using the unfitted scheme one is now 
much greater than either scheme two or three. A third such comparison 
with = 0.4 is shown in Figs. 9.19 and 9.20. Figures 9.19b and 9.19c 
show Mach contours generated from the results of scheme two on the 
coarse and fine grids respectively. The contours shown in Fig. 9.19a 
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using scheme one and Fig. 9.19c -using scheme two on the fine grid both 
show a supersonic zone extending across the channel just downstream of 
the expansion corner which suggests that the flow is chocked in the 
channel and then expands supersonically as in a converging-diverging 
nozzle. The comparison of "M" vs. £ in Fig. 9.20 shows that the flux 
vector splitting results show the least separation whereas the results 
from scheme one are the most separated. 

The entropy produced by flow through the shock is shown in Figs. 

9.21a-9.21c. These results were obtained using scheme two with the 

upper boundary = 20. In Fig. 9.21a, the lower boundary Is at 2.0 

and the maximum value of the entropy is 0.0878 just downstream from the 

shock along the wall. In Fig. 9.21b, the lower wall n 0 is at 1.6 and 
the maximum value of S is 0.126 and the entropy gradients are somewhat 
larger than in Fig. 9.21a. In Fig. 9.21c, the lower wall n 0 is at 1.2 
and the maximum value of S is 0.145. It can be seen that there is a 
recirculation region and the entropy gradients are again greater than in 
the previous two cases with most of the entropy being produced very 
close to the lower boundary. These results correlate with Fig. 9.10, 
which shows the vorticity level increasing as n 0 is reduced and with 
Crocco's theorem which states the vorticity will exits in a flow if 
gradients of entropy and/or enthalpy also exist [60], 

Figures 9.22 through 9.25 show the results of the flux splitting 
scheme three with the upper boundary at 40 and the lower boundary being 
reduced from 0.4 to 0.2 to 0.1. In Figs. 9.22a and 9.22b, the Mach 
contours and sonic line are shown for n 0 s 0.4 and the large separated 
zone downstream of the step and the shock at the expansion corner are 
The Mach contours for n Q « 0.2 in Fig. 9.23a and the 


evident. 
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corresponding sonic line plot “in Fig. 9.23b show a much smaller 
supersonic zone at the corner and a slightly larger supersonic zone 
which is apparently not terminated by a shock just slightly downstream 
and above it. In Fig. 9.24a, the Mach contours for n Q = 0*1 are shown 
and the corresponding sonic zones are given in Fig. 9.24b. It can be 
seen that the supersonic zone at the corner has nearly disappeared and 
the supersonic bubble out in the flow field has grown. The entropy 
contours shown in Fig. 9.24c and the stream function plot in Fig. 9.24d 
show clearly the separation region with the entropy now associated with 
the rotating vortex. Figure 9.25 shows the portion of the grid around 
the expansion corner and the resolution which is achieved at the corner 
by the special clustering. 

A similar procedure was done using the same scheme but with the 
upper boundary instead at * 20. The results at n Q = 0.1 are shown 
in Figs. 9.26 and 9.27. Figure 9.26a shows the Mach contours and it is 
apparent that a shock has formed which extends across the channel 
downstream of the expansion corner as was noted earlier. There is a 
large supersonic zone as is shown in Fig. 9.26b, which is terminated by 
the shock. The flow has separated from the corner and there is entropy 
trapped inside the recirculation region as shown in Fig. 9.26c. Figure 
9.27 shows the Mach numbers multiplied by the sign of the contravarient 
velocity in the streamwise direction. The location of the shock at the 
expansion corner and at the top wall are clearly evident. It is also 
apparent that the shock at the expansion corner is still fairly strong 
and may not be going to zero as n 0 approaches 0. 









141 




143 



Fig. 9.25 Grid Detail at Corner. 
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Fig. 9.26(c) Entropy Contours 
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9.1.2 Viscous Results — 

The Navier-Stokes scheme four has run on the rearward facing step 
under a variety of conditions. The purpose of the runs was to determine 
how the viscous results might differ from the inviscid results. No 
attempt was made to duplicate all of the various positions of the lower 
wall n Q and the upper wall due to the time and computational expense 
that would be necessary to do this and also because, as will be seen, 
the results that were obtained were so different from the inviscid 
results that this was felt by the writer to be unnecessary. 

The first set of runs was done using a Reynolds number of 10,000 
with the upper boundary at nj = 40 and the lower boundary at n 0 = 2. A 
grid size of 81x61 was used and the first point off the wall was at 
n * 0.01. A constant global time step of 0.2 was used at the solution 
was run for 3000 iterations and plots obtained at 1500, 2000, 2500 and 
3000 iterations. Figure 9.28 shows the results at 1500 Iterations. 
From the Mach contours in Fig. 9.28a, it is evident that the boundary 
layer which begins upstream of the corner separates at the corner and 
that there is no shock there. The solution throughout the flow field is 
fairly smooth and the freestream upper boundary condition does not 
appear to be causing any flow abberations. Figure 9.28b shows the 
pressure contours and the presence of at least one vortex is apparent 
from the circular pattern with a low pressure center. Figure 9.28c 
shows the velocity vectors and this large vortex is apparent along with 
two other smaller vorticies upstream. The wall shear stress plot in 
Fig. 9.28d provides further evidence of flow reversal along the lower 
wall between approximately £ = -4.6 and 5 = -1 and also in the vicinity 
of 10 <£< 17. 




(c) Velocity 


(d) 


Wall 


100 


50 


-50 


-100 


Shea 


Fig. 9.28 (co 


WSS 


149 


The results of 2000 iterations are shown in Figs. 9.29a through 

9.29d. A comparison between Figs. 9.28a and 9.29a shows that the flow 
field has changed greatly in that the extent of the separation zone has 
grown and there now appears to be more disruption of the flow field 
above the separated region. The pressure contours in Fig. 9.29b an the 
velocity vector plot in Fig. 9.29c show the presence of two large 

vorticies downstream of the corner. The wall shear stress plot in Fig. 
9.29d shows a third reverse flow region in the region of -5<5<-l ahead 
of the step. 

Ther results at 2500 iterations are shown in Figs. 9.30a through 

9.30d and it is apparent that the location of the vortices has again 

changed which indicates that either the flow has not reached a final 
steady state or the flow may be unsteady. The writer considers the 
later possibility to be the more likely. One large vortex is present 
just downstream from the corner and the writer believes it is possible 
that this may be the small vortex which was just beginning to form at 
2000 iterations and has now moved further downstream and grown. At the 
same time, the two large vortices that were present at 2000 iterations 
have been swept out of the computational region. Again, from the wall 
shear stress plot in Fig. 9.30d it is apparent that a smaller reverse 
flow region appears to exist inside the boundary layer ahead of the 
step. This region has not changed location significantly from the 
previous two iterations and it may be a permanent feature that has 
established itself and may be giving rise periodically to shed vortices 
which are then swept downstream. 

The results at 3000 iterations are shown in Figs. 9.31a through 
9 . 3 Id and again the flow field appears to have changed. Although the 
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Mach contours in Fig. 9.31a do resemble the Mach contours in Fig. 9.30a, 
the pressure contours in Fig. 9.31b now show a fairly strong adverse 
pressure gradient ahead of the corner which was beginning to become 
apparent in Fig. 9.30b. Figure 9.31b shows the presence of four 

circular shaped pressure regions of alternating high and low pressure. 
The velocity vector plots in Fig. 9.31c show two clockwise rotating 
vorticies which appear to correspond in location to the two low pressure 
regions in Fig. 9.31b. The high pressure regions appear to correspond 
to the regions where the flow velocities are low just downstream of the 
step and in between the two vorticies. The wall shear stress plot in 
Fig. 9.31d again shows a small recirculation zone ahead of the step 
which must lie inside of the boundary layer. 

Two runs were also done at a Reynolds number of 100,000 .and the 
results are compared in Figs. 9.32 and 9.33. The first run was done 
using an 81x61 grid with a freestream Mach number of 0.865 using a 
nondimensional time step of 0.02. The flow was assumed to be laminar 
and the turbulent terms were turned off. The Mach contour, pressure 
contour, wall shear stress, and velocity vectors after 2000 iterations 
are shown in Figs. 9.32a-9.32d respectively. It can be seen that the 
results are not appreciably different than the results at the lower 
Reynolds number shown previously. Ther are numerous recirculation 
vorticies present in the flow as can be seen from the velocity vector 
plot and there are no shocks in the flow. The second set of results was 
obtained using the same Reynolds number and time step but a 161x61 grid 
and the flow was assumed to be turbulent. The results after 2000 

iterations are shown in Figs. 9.33a-33d. From the Mach and pressure 
contours and the wall shear stress plot, it is apparent that a strong 
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shock now exists across the channel similar to some of the inviscid 
cases. Furthermore, the velocity vector plots show only a single 
recirculation vortex which is also similar to the inviscid cases. It 
thus appears that putting the turbulence model into use results in a 
flow pattern that strongly resembles some of the inviscid results. 

9.2 Problem Two 


9.2.1 Inviscid Results 

The second configuration, a "bump" inside a channel, was treated in 
a manner similar to the rearward facing step. The three inviscid 
schemes were all run on 81x41 grids and the shock fitting scheme two was 
also run on a fine grid 161x81 for comparison. The Navier-Stokes code 
was run on an 81x61 grid and the results from using this code is 
discussed in the next section. In all cases, the upper wall nj was set 
at 12.0 and the left and right boundaries were set at £ = -10.0 and 
+10.0 respectively. The lower wall n 0 was varied from a maximum of 1.0 
to a minimum of 0.2 in the inviscid cases. The back pressure boundary 
condition was that the nondimensional pressure p 1 is 0.5. 

Figure 9.34 shows the Mach number on the upstream side of the shock 
foot for the various cases which were run. The shock capturing results 
from scheme one again show that as the lower wall n Q is reduced and the 
expansion becomes greater that the Mach number ahead of the shock 
increases much faster than the other schemes which is a result of the 
nonconservative nature of the scheme. The shock fitting results on both 
the coarse and fine grids agree fairly closely with the Mach number 
increasing to about 2.2 as the lower wall n 0 is reduced. The flux 
vector splitting scheme three results show lower Mach numbers than the 
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two previous schemes and this is-consistent with the results of problem 
one. The maximum Mach number achieved with scheme three is 1.92 at n 

o 

equal to 0.6. 

Figure 9.35 shows the shock locations for the four sets of results. 
The scheme one shock location is consistently greater than 1.4 and 
reaches a maximum of 1.61 at n 0 = 0.2. It should be noted that since 
the grid is fixed in cases one and three, the shock location shifts of 

increment of one grid point and cannot achieve intermediary values. The 

shock fitting scheme two results show that the shock location starts out 
at the same location as scheme one but that as the lower wall n Q is 
reduced, the s is gradually reduced to about 0.7 at a lower wall n 0 of 
0.2. The results of scheme three are similar at n Q equal to 0.2, the 
S is 0.53. 

Figures 9.36 and 9.37 show the entropy and vorticity downstream of 

the shock adjacent to the wall. In Fig. 9.36, it can be seen that as 

the lower wall is reduced the entropy produced by flow through the shock 
tends to increase in all cases. The shock captured solutions of scheme 
one show less entropy than the results of the shock fitting scheme or 
the flux vector splitting scheme three. The results of schemes two on 
both the coarse and fine grids and scheme three agree fairly closely. 
In Fig. 9.37, it can be seen that in all cases the vorticity is 
increased as the lower wall is reduced. However, the actual levels 
of vorticity vary depending on the scheme. The flux vector splitting 
scheme three shows the least the amount of vorticity present in the 
flow; however, this may be due to the fact that the function values in 
scheme three are cell volume average instead of point values and this is 
reflected in the way the values are reported out. 
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The minimum Mach numbers multiplied by the sign of the 
contravarient velocity component along the lower wall for the various 
cases are shown in Fig. 9.38. The results all indicate that as n Q is 
reduced, the value of M decreases and that the shock fitting results on 
the coarse grid were the first to separate which is not the same as 
problem one when the scheme one results were the first to separate. In 
general, however, the results in Fig. 9.38 show a fairly high degree of 
consistency. The 5 location of the flow reattachment shown in Fig. 9.39 
indicate that the results of scheme one generally show a smaller 
separation zone than either shock fitting scheme two or the implicit 
shock fitting scheme three. 

A comparison of the results of the three schemes for the case with 
the lower wall n 0 equal to 1.0 is shown in Figs. 9.40a-9.40d. From the 
Mach number contours which are shown, it can be seen that the results 
are all similar. The fine grid shock fitting results are not as well 
converged as the coarse grid results since both runs were for 1000 
iterations and the fine grid case would need probably 2000 iterations 
(which would take 8 times as long since there are roughly four times as 
many points as the coarse grid) to get to the same level of convergence. 

A similar comparison is shown in Figs. 9.41-9.41d. The Mach 
contours produced by scheme one are very similar to those produced by 
scheme two and, in fact, if the two plots are overlaid, the contours are 
seen to be nearly identical throughout the flow field except close to 
the shock. The flux vector splitting Mach contours show a smaller sonic 
region but overall are similar to the results of schemes one and two. A 
comparison of the Mach number multiplied by the sign of the contra- 
varient U shows that the jump through the shock is significantly less 
with scheme three than with schemes one and two. 
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The next set of results which are shown In Figs. 9.42-9.43 compare 

the coarse grid shock fitting results to the flux vector splitting 

results for the case when n is 0.2 which is the most severe case which 

o 

was run. The Mach, entropy, and stream function plots for the shock 
fitting results are shown in Figs. 9.42a-9.42d respectively, and the 
same plots for the flux vector splitting results are shown in Figs. 

9.43a-9.43d respectively. A comparison of the Mach contours show that 
the flux vector splitting results shows a smaller sonic region and that, 
in general, the contours are smoother. A comparison of the entropy 
contours shows that the shock fitting scheme is particularly good at 
preserving the correct entropy jump through the shock with zero entropy 
in the flow which has not been passed through the shock. The vorticlty 
contours show that the shock fitting scheme results in some vorticity 
being present in the flowfield above the shock where it should not be 
and this is possibly due to lack of grid smoothness in this region. The 
vorticity plot Fig. 9.43d shows high gradients of vorticity near the 

lower boundary which may be due to some error in the way these values 
were computed although this has been checked and is felt not to be the 
case. The comparison of the Mach number along the lower boundary 
multiplied by the contraviant U in Figs. 9.42b and 9.43b shows again 

that the flux split results have a smaller jump in M through the shock 

than the shock fitting results. 

9.2.2 Viscous Results 

The Navier-Stokes code was run for 2000 iteration at a Reynolds 
number of 10,000 and the results are shown in Figs. 9.44 and 9.45. The 
flow was assumed to be laminar and calculated without the turbulence 
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model. In Fig. 9.44a, the Mach and pressure contours after 1500 
iterations are presented. The lower wall n Q is equal to 1.0 and it can 
be seen that the flow has separated downstream of the step with at least 
one vortex present. The results at 2000 iterations in Fig. 9.45 also 
show the flow to be separated but the Mach and pressure contours are 
different from those at 1500 iterations indicating that the flow is 
unsteady. An examination of the convergence history reveals that after 
1500 iterations, the computations do not converge any further, which is 
a further indication that the flow is unsteady. The viscous results are 
thus very different from the inviscid results and show that the inviscid 
separation phenomenon is not a good predictor as to when viscous 
separation will occur and that viscous calculations are needed. 

9.3 Problem Three 

The problem of flow past a circular arc (occasionally referred to 
as "Ni's bump") was solved using the shock capturing scheme 1 and the 
flux vector splitting scheme 3. An effort to do shock fitting was made 
using two different approaches but both yielded unsatisfactory results. 
The reasons for this are discussed in this section and a possible 
solution is described. The difficulties are the result of the fact that 
unlike the previous two configurations, this configuration has sharp 
stagnation corners where the flow velocity goes to zero. It is, 
therefore, desirable to cluster grid points near the corners in order to 
capture the high flow field gradients that exist in these regions. 
Since the shock fitting, which was attempted, moves the grid as the 
shock moves, special attention must be given to insure that the shock 
movement does not move the clustered grid away from these corners. 
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Figures 9.46a shows the first grid which was used to solve the 
problem using scheme one. The grid was developed using bilinear 
interpolation and the elliptic smoothing described in Chap. 2. The grid 
achieves good clustering at the stagnation corners and over the bump 
through the third order polynomial distribution which was used. The 
nondimensional back pressure p/p 0 for this case was 0.74. The scheme 
was run for 1000 iterations and the resulting Mach contours are shown in 
Fig. 9.46b. It can be seen that the flow is transonic with a shock 
present over the bump. Unfortunately, there are some wiggles present in 
the flow field emanating from the corners which are apparently due to 
the grid clustering which leads to high gradients in the metric terms. 

The writer was not able to do the shock fitting on this grid since 
it was generated using algebraic and elliptic methods instead of 
conformal mapping and scheme two relies on the existence of a conformal 
mapping to get the new transformation metrics and x^ each time the 
grid was moved. In principal, it is possible to do this by interpolat- 
ing the metrics from the original grid (since they are known) to the new 
shock fitted grid once the shock position is known and the writer tried 
this approach but it was not successful. 

The writer then learned that a conformal transformation exists for 
this configuration, therefore, it ought to be possible to use the same 
shock fitting program on this case. The transformation was then used to 
develop the 81x41 grid shown in Fig. 9.47a, which can be seen clusters 
the grid points over the bump but not in the corners. Scheme one was 
then run for 2000 iterations on this case and the resulting Mach 
contours are shown in Fig. 9.47b. Although the back pressure is the 
same as in the previous run, the shock on top of the bump is much 
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weaker. A plot of the Mach numbe-p- along the lower wall is shown in Fig. 
9.47c and it can be seen that there is a fairly weak shock at approxi- 
mately 5 = 1.0. It is also evident that true stagnation conditions are 
not achieved at the corners of the bump as they should be. It was 
throught that a higher grid resolution might alleviate these problems so 
the writer then ran the same case using a 161x41 grid. The resulting 
Mach contours are shown in Fig. 9.48a and the Mach number along the wall 
is shown in Fig. 9.48b. It is interesting to note that the solution 
near the lower boundary has high gradients in the Mach numbers and that 
the velocities near the lower boundary are lower than in the flowfield 
just off the wall. This is believed to be due to the fact that there is 
a discontinuity in the metric terms at the corners wich results from the 
way the grid has been generated. 

In order to obtain a shock captured solution which might facilitate 
shock fitting better than the previous case, it was decided to move the 
upper boundary further away and to increase the back pressure to a 
nondimensional p/p 0 of 0.54. Scheme one was run on a 81x41 grid and the 
resulting Mach contours are shown in Fig. 9.49a and the vorticity 
contours in Fig. 9.49b. The Mach contours show a strong shock on the 
leeward side of the bump. The vorticity contours in Fig. 9.49b show 
high vorticity along the shock and also at the corner on the upstream 
side of the bump. The vorticity upstream of the shock is not physically 
correct and indicates again that the scheme does not do well in such 
regions. The shock fitting scheme two was tried on this case but with 
little success. The Shock would become unstable after the code had run 
for less than 100 iterations and the code would encounter execution 


errors and stop 
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The grid was next refined to a 161x41 grid and the shock capturing 
scheme one tried again. The resulting Mach and entropy contours are 
shown in Figs. 9.50a and 9.50b respectively. It was found that refining 
the grid again resulted in a poor solution near the lower boundary due 
to the dicontinuity in the grid metrics at the corners. The entropy 

contours in Fig. 9.50b show that there is an entropy layer produced 

which starts at the front corner. This is an unphysical result and the 
poor quality of the solution made it not feasible to go on to the shock 
fitting. 

The flux vector splitting scheme three was run on the 81x41 grid 

and the resulting Mach contours are shown in Fig. 9.51. The Mach 

contours appear smooth and there is no problem at the boundaries. The 
writer feels that this is a much better solution and although the shock 
is not resolved as well as it might be, the case using shock fitting; 
the solution was obtained without any special treatment which would have 
been necessary to get a shock fitted solution. 

9.4 Problem Four 

The grid for the NACA 0012 airfoil was produced by the elliptic 
procedure described in Chap. 2. The grid covered only the top half of 

the airfoil and used 91x33 points in the e; and n coordinated directions. 

\ 

The outer boundary was placed ten chords away from the airfoil surface 
and the downstream boundary was ten chords lengths away. There are 73 
points along the top surface of the airfoil. The grid is shown in Fig. 
9.52 in three perspectives ranging from far away to close to the leading 
edge. It can be seen that the resolution near the leading and trailing 
edges is reasonably dense and that the orthogonality of the grid at the 
airfoil surface is good. 
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(a) 91x31 C-Grid 



(b) Detail over airfoil 
Fig. 9.52 Grid for NACA 0012 Airfoil. 
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Scheme one was specially modified to run this particular case and 
the freestream Mach number used was 0.8, as this is a good transonic 
test case. The scheme was allowed to run for 3000 iterations which 
resulted a drop in the residual of about 2.5 order of magnitude. The 
Mach contours are shwon in Fig. 9.53. Although the solution may appear 
to be reasonable, there are again kinks in the contours near the 
trailing edge which may be due to lack of grid smoothness. Also, the 
shock location and strength are not correct. Steger has computed this 
case and gives the shock location in terms of x/L where L is the chord 
as being 0.6 and the Mach number just before the shock as 1. The 
results which were computed using scheme one show the shock at 0.73 of 
the chord and the Mach number as 1.405. The solution is, therefore, not 
correct and either shock fitting or the flux vector splitting code is 
necessary for better results. 

The writer modified scheme two to do this case, but it was not 
stable and the grid developed extreme skewness as a result use of a 
fifth order polynomial function to distribute the grid points with the 
shock moving. The writer feels that a shock fitting scheme which has 
the grid stationary and simply extrapolates to the shock from the flow 
field would work much better than the existing shock fitting scheme 
which attempts to adapt the grid to the moving shock. Such a scheme is 
being developed by the writer but at the time of this writing, it has 
not yet been completed or tested. 




Chapter 10 


CONCLUSIONS 

This study was originally initiated as an investigation into the 
inviscid separation phenomenon. Since that time, there have been 
several studies on this topic published and, as yet, reseachers do not 
seem to agree on the significance of these inviscid results. The 
questions concerning inviscid separation appear to belong to three 
separate but related categories which can be stated as: 

1. When are numerically computed solutions which exhibit flow 

separation valid solutions to the governing equations and are they 
unique? 

2. In the case of computed solutions where closed recirculation regions 
are present, what is the source of the vorticity inside the eddies? 

3. What is the physical relevance of an inviscid separated flow to the 
viscous case? 

The first question has been investigated by various researchers 
either by grid refinement procedures or by careful control of parameters 
such as the amount of artificial dissipation which is added to a scheme 
to provide stability [15, 17, 18, 20, 22]. In general, it appears that 
while in some cases the computed inviscid solutions are valid and 
reasonable, it is also possible to pcoduce inviscid separation if 
excessive artificial dissipation is present or if very coarse grids are 
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used. The question of uniqueness has not been adequately examined, 
although it is known that in some cases nonunique solutions to the Euler 
equations do occur. Since numerical solutions are by nature approximate 
solutions and there exists few rigorous proofs that the solutions 
obtained by any given numerical method are in fact the correct 
approximations to the problem (Including boundary conditions) being 
solved, it is possible that this issue will never be completely 
resolved. 

There does not, as yet, appear to be any clear answer to the second 
question. The most common explanation for the presence of vorticity is 
that flow through curved shock waves results in entropy gradients and 
thus by Crocco's theorem, vorticity. Also, since in nearly all cases 
reported, the inviscid separated zone are "captured" during the 
iterations procedure by the numerical algorithm and the initial assumed 
flow field is free of such separated zones, it can be argued that the 
unsteady transient flow leads to entropy and vorticity in the part of 
the flow field that is eventually "pinched-off " and trapped inside the 
eddy. It should also be noted that while in theory, there exists a 
dividing stream line between the trapped eddy and the outside flow field 
and once formed, vorticity and entropy cannot enter the eddy, in actual 
numerical computations the flow inside the eddy is coupled to the flow 
outside the artificial viscosity and the boundary between the two zones 
(which is smeared) may allow the total amount of entropy and vorticity 
inside the change. 

The third question may be the most important of all, since usually 
the purpose of doing an inviscid computation is to provide some insight 
as to how an actual viscous flow would behave without the expense of 
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actually doing viscous calculations. Since real fluids contain 
viscosity, the first two questions, while of theoretical interest, do 
not pertain to the numerical simulation of viscous flows. 

This investigation had many different aspects to it, but the 
primary effort involved an attempt to solve the transonic flow past 
various two-dimensional configurations using various numerical 
procedures. The four configurations which were chosen included a 
rearward facing step, a "bump" in a channel, a circular arc airfoil, and 
a NACA 0012 airfoil. For the first two configurations, the position of 
the lower boundary was varied to study the effect of flow separation a 
the curvature in the expansion region increased. The last two 
configurations were used to determine how the numerical schemes could 
handle somewhat more complicated geometries. 

The grids for these configurations were produced using either 
conformal mapping or algebraic and elliptic procedures. The conformal 
mapping procedure has the advantage that it is fast and straightforward 
and produces grid which are perfectly orthogonal as well as satisfy the 
Cauchy-Riemann equations. Also when doing shock fitting, it is 
necessary to recompute the metric terms each time the grid is moved and 
this is done very easily if there exists a simple transformation which 
can be used. The disadvantage is that a different mapping must be found 
for each configuration which greatly restricts the geometrical 
versatility of the method. Finally, the method cannot be extended to 
three-dimensional regions in the general case. The algebraic and 
elliptic methods were used for the last two problems and worked well. 
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The numerical schemes which-were used included three Euler solvers 
and a Navier-Stokes solver. The three Euler solvers included an 
explicit Lambda scheme, a shock fitting version of this scheme and a 
flux-vector splitting scheme. The first two schemes were coded by the 
writer and the shock fitting version included some features developed by 
the writer to enable the code to do imbedded shock fitting for the 
regions which were considered. These included the use of polynomial 
functions to smoothly blend the end of the shock to the grid line which 
extended from the end of the shock to the upper boundary and the use of 
least squares fitting to smooth the shock acceleration, velocity and 
position as it was being computed in order to overcome stability 
problems which were encountered. Even so, the method was not successful 
in all cases. The first two schemes were also vectorized for the NASA 
Langley VPS 32 supercomputer which greatly increases their execution 
speed. The flux vector splitting scheme proved to be the most robust 
and was easy to use. Generally, the convergence rate was also higher 
with this code than with the explicit schemes although since it was not 
vectorized, the overall execution time to get a converged solution was 
of the same order of magnitude as the explicit schemes. The Navier- 
Stokes code was also robust and was vectorized for the VPS 32. 

The three inviscid schemes were run on configurations one and two 
under a variety of conditions which were intended to examine the effects 
of changing the position of the lower boundary (on configuration one the 
upper boundary was also tested in three positions). The Navier-Stokes 
was run under more limited sets of conditions since the writer felt that 
the runs which were conducted were sufficient for comparison purposes. 


197 


Schemes one and three were nin on configuration four and an attempt 
was made to do the shock fitting on this configuration but this was 

unsuccessful. A new way to do the shock fitting, which does not move 
the grid with the shock, is currently being developed but is not 
available for use as of this writing. This method was not felt to be 
necessary when this study was initiated and therefore was not included 
in the proposal for this research. Scheme one was run on the NACA 0012 
configuration four, but the position and strength of the computed shock 
were not correct. Shock fitting has also been tried on this case but 
again was not successful. However, the flux vector splitting scheme has 

been run on this case extensively by von Lavante and the results as 

reported by this scheme agree well with other reported results. 

The major conclusions, which the writer believes can be made from 
this investigation, are summarized below. 

1. Conformal grid generation worked well for the simple two- 

dimensional configurations which were generated as test cases. It 
is fast and produced perfectly orthogonal grids. 

2. For configurations three and four, algebraic and elliptic grid 

generation was used successfully. This method has the advantage 
that a grid can be produced for virtually any configuration. The 

algebraic method is fairly fast but must be tuned to get good 

results. The elliptic method will produce good orthogonality at 
the boundaries provided the source terms are carefully 

controlled. The elliptic method requires the solution of a coupled 
set of POE's and, therefore, results much more computer time than 
either the algebraic or conformal methods. 
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3. The shock capturing Lambda .scheme is an accurate scheme which was 
easy to code and is based on characteristic theory. The scheme 
gave results which were in good qualitative agreement with the 
other inviscid schemes used. Unfortunately, since the equations 
are solved in nonconservative form, the results at shocks were 
generally not correct. Also, the scheme converged very slowly due 
to the explicit CFL limitation. 

4. The shock fitting scheme developed by the writer overcame the 
disadvantage of the previous scheme by explicitly solving the 
Rankine-Hugoniot jump relations through the shock. An adaptive 
grid which moved with the shock was used to obtain a high 
resolution of the correct shock shape and position. Vectorization 
of the scheme was fairly straightforward and made the computations 
fairly efficient despite the restrictive CFL limitation. The major 
disadvantage is that the scheme (at least as it was used by the 
writer) requires a conformal transformation to exist to enable a 
rapid evaluation of the metric terms which must be updated each 
time the grid is moved. For general configurations where no such 
transformation exists, shock fitting is still possible but would be 
more difficult to implement. 

5. The flux vector splitting scheme was a robust code and achieved the 
correct jump relations through the captured shocks without the use 
of fitting. The scheme converged faster than the explicit codes 
since the code is implicit which allows a much higher CFL number to 


be used. 
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6. Two configurations, a rearward facing step and a bump inside a 
channel, were studied for flow separation. The inviscid cases were 
run using various positions of the lower boundary to compare when 
separation occurred. It was found that when the lower boundaries 
were such that the channel transition were very gradual, no 
separation occurred. However, as n 0 was brought closer to zero, 
flow separation would eventually occur and it appeared to be 
related to entropy gradients in the flow field downstream of the 
shock which produced vorticity in the flow. 

7. Grid refinement was done using the vectorized shock-fitting scheme 
and the results did not indicate that the flow was less likely to 
separate than on the coarse 81x41 grid. 

8. As n Q is reduced to zero in the rearward facing step, it appears 
form the results using the flux-vector splitting scheme that the 
shock at the corner in the inviscid cases eventually is reduced to 
a negligible size and the flow simply separates smoothly from the 
corner with a single recirculation eddy. This is in agreement with 
previous results obtained by Jameson [61]. This was not, however, 
observed with the shock fitting code as it was not possible to get 
results with lower than 0.4 since the resolution of the shock- 
fitting grid was not adequate at the corner. 

9. The viscous results which were obtained were significantly 
different from the inviscid results in that flow separation 
occurred due to the effect of the adverse pressure gradients on the 
boundary layer and the flow was steady. In the inviscid cases, 
separation appears to be due to vorticity at the wall produced by 
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flow through the shock {wWch was not present In the viscous 
cases). The viscous results generally showed an unsteady shedding 
of vorticities from the leeside of the corner or the bump with 
multiple recirculation zones; whereas, the inviscid results showed 
no separation for the same position of the lower boundary. When 

the lower boundary n 0 was reduced, a single recirculation zone 
formed in the inviscid cases which was steady and was not shed. 

10. Two additional configurations were investigated with the shock 
capturing schemes and an attempt was made to do shock fitting. The 
first was a simple circular arc airfoil and the second was a NACA 
0012 airfoil. The shock fitting scheme was attempted on the 
circular air airfoil but the conformal grid which was generated had 
a discontinuity In the metric terms at the corner of the arc with 
the lower boundaries which resulted in the production of an entropy 
layer emanating from the front stagnation corner which made the 
shock fitting Lambda scheme unstable. The flux vector splitting 
scheme was run successfully on this configuration. 

11. The NACA 0012 case was run with the shock capturing scheme one but 
the shock position and strength were nor correct and shock fitting 
was tried but so far has not been successful. This is due in part 
to the need to adjust the grid by a fifth order polynomial to keep 
adequate resolution at the leading and trailing edges. A different 
type of shock fitting which holds the grid stationary while the 
shock is moved may be necessary to avoid this problem. 
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